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Advective  diffusion  of  passive  tracers  by  fluid  flow  plays  a  key  role  in  the  transport 
of  salt,  heat,  buoys  and  markers  in  geophysical  flows,  in  the  dispersion  of  pollutants 
and  trace  gases  in  the  atmosphere,  and  even  in  the  motion  of  sea  ice  floes  influenced 
by  winds  and  ocean  currents.  The  long  time,  large  scale  behavior  of  such  systems  is 
equivalent  to  an  enhanced  diffusion  process  with  an  effective  diffusivity  tensor  D*. 
Over  twenty  five  years  ago,  a  Stieltjes  integral  representation  for  the  homogenized 
tensor  was  derived,  involving  the  spectral  measure  of  a  self-adjoint  operator  that 
depends  on  the  flow  characteristics.  However,  analytical  calculations  of  D*  have 
been  obtained  for  only  a  few  simple  flows,  and  numerical  computations  of  the  effective 
behavior  based  on  this  powerful  representation  have  apparently  not  been  attempted. 
We  overcome  these  limitations  by  providing  a  mathematical  foundation  for  rigorous 
computation  of  D*.  We  explore  two  different  approaches  and  for  each  approach  we 
derive  Stieltjes  integral  representations  for  the  symmetric  and  antisymmetric  parts 
of  D*,  involving  a  spectral  measure  of  a  self-adjoint  random  operator.  In  discrete 
formulations  of  each  approach,  we  express  the  spectral  measure  explicitly  in  terms 
of  standard  (or  generalized)  eigenvalues  and  eigenvectors  of  Hermitian  matrices.  We 
develop  and  implement  a  numerically  efficient  projection  method  that  significantly 
reduces  the  dimension  of  the  random  matrix  in  computations  of  spectral  statistics. 
We  use  this  method  to  compute  the  effective  behavior  for  several  model  flows  and 
relate  spectral  characteristics  to  flow  geometry  and  transport  properties. 


2 


*  University  of  Utah,  Department  of  Mathematics,  155  S  1400  E,  RM  233,  Salt  Lake  City,  UT  84112-0090, 


USA 


3 


I.  INTRODUCTION 

The  enhancement  in  diffusive  transport  of  passive  scalars  by  complex  fluid  flow  plays  a 
key  role  in  many  important  processes  in  the  global  climate  system  [55]  and  Earth’s  ecosys¬ 
tems  [14] .  Advection  of  geophysical  fluids  intensifies  the  dispersion  and  large  scale  transport 
of  heat  [36],  pollutants  [7,  12,  48],  and  nutrients  [14,  21]  diffusing  in  their  environment. 
Advective  processes  also  enhance  the  large  scale  transport  of  plankton  [21],  which  is  an 
important  component  of  the  food  web  that  sustains  life  in  the  polar  oceans.  The  transport 
of  vast  sea  ice  floes  in  the  polar  oceans  is  driven  by  a  seasonally  and  regionally  changing 
balance  in  oceanic  and  atmospheric  forces  [28,  55].  These  forces  can  enhance  the  transport 
of  sea  ice  floes  by  eddie  fluxes  in  the  atmosphere  and  ocean  currents  [28,  44,  45,  56]. 

Complex  interactions  between  shearing  ocean  waves,  tidal  currents,  and  wind  drift,  for 
example,  gives  rise  to  complex,  behavior  of  the  flow  fields  [11,  27,  59].  It  was  discovered  in 
the  early  1900s  [53]  that  complex  fluid  flows  transport  passive  scalars  in  much  the  same  way 
as  molecular  diffusion.  The  mathematical  description  of  this  phenomenon  [54]  demonstrated 
that  the  long  time,  large  scale  behavior  of  a  diffusing  particle  or  tracer  being  advected  by  an 
incompressible  fluid  velocity  field  is  equivalent  to  an  enhanced  diffusion  process  with  an  effec¬ 
tive  diffusivity  tensor  D*.  Describing  the  enhancement  of  the  effective  transport  properties 
by  fluid  advection  is  a  challenging  problem  with  theoretical  and  practical  importance  in  many 
fields  of  science  and  engineering,  ranging  from  turbulent  combustion  [1,  10,  43,  52,  57,  58] 
to  mass,  heat,  and  salt  transport  in  geophysical  flows  [36]. 

A  broad  range  of  mathematical  techniques  have  been  developed  that  reduce  the  analysis 
of  complex  fluid  flows,  with  rapidly  varying  structures  in  space  and  time,  to  solving  averaged 
or  homogenized  equations  that  do  not  have  rapidly  varying  data,  and  involve  an  effective 
parameter  [9,  15,  16,  29-31,  41,  58].  Homogenization  of  the  advection-diffusion  equation  for 
passive  scalar  transport  by  random,  time-independent,  mean-zero  fluid  velocity  fields  was 
treated  in  [31].  This  reduced  the  analysis  of  advective  diffusion  to  solving  a  diffusion  equation 
involving  a  homogenized  temperature  and  a  (constant)  effective  diffusivity  tensor  D*. 

An  important  consequence  of  this  analysis  is  that  the  effective  diffusivity  D*  is  given  in 
terms  of  a  random  “cell  problem”  involving  a  curl-free  stationary  stochastic  process  [31], 
which  satisfies  a  steady  state  diffusion  equation  involving  a  skew-symmetric  random  matrix 
H  [2,  3,  15,  16].  By  adapting  the  analytic  continuation  method  of  homogenization  theory 
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for  composite  materials  [19],  it  was  shown  [2,  3]  that  the  cell  problem  could  be  written  as  a 
resolvent  formula  involving  a  self-adjoint  random  operator  acting  on  the  Hilbert  space  of  curl- 
free  vector  fields.  This,  in  turn,  led  to  a  Stieltjes  integral  representation  for  the  symmetric 
part  of  D*,  involving  the  Peclet  number  Pe  of  the  flow  and  a  spectral  measure  of  the  operator. 
A  key  feature  of  the  integral  representation  for  D*  is  that  parameter  information  in  Pe  is 
separated  from  the  geometry  of  the  fluid  velocity  held,  which  is  encoded  in  the  spectral 
measure  through  its  moments.  This  parameter  separation  has  led  [2,  3]  to  rigorous  forward 
bounds  for  the  diagonal  components  of  D*. 

The  mathematical  framework  developed  in  [31]  was  adapted  [41]  to  the  case  of  a  periodic, 
time-dependent  fluid  velocity  held  with  non-zero  mean.  It  was  shown  [41]  that,  depending 
on  the  strength  of  the  fluctuations  relative  to  the  mean  how,  the  effective  diffusivity  tensor 
D*  can  be  constant  or  a  function  of  both  space  and  time.  When  D*  is  constant,  only  its 
symmetric  part  appears  in  the  homogenized  equation  as  an  enhancement  in  the  diffusivity. 
However,  when  D*  is  a  function  of  space  and  time,  its  antisymmetric  part  also  plays  a 
key  role  in  the  homogenized  equation.  Based  on  [8],  the  cell  problem  associated  with  a 
time-independent  how  was  transformed  [41]  into  a  resolvent  formula  involving  a  self-adjoint 
operator,  acting  on  a  Sobolev  space  [17,  33]  of  spatially  periodic  scalar  fields,  which  is  also 
a  Hilbert  space.  This,  in  turn,  led  to  a  discrete  Stieltjes  integral  representation  for  the 
antisymmetric  part  of  D*,  involving  the  Peclet  number  of  the  how  and  a  spectral  measure 
of  the  operator. 

Such  methods  have  been  extended  to  steady  hows  where  particles  diffuse  according  to 
linear  collisions  [42],  solute  transport  in  porous  media  [8],  anelastic  (weakly  compressible) 
hows  [32],  and  to  the  setting  of  a  time-dependent  fluid  velocity  held  [4,  38].  All  these 
approaches  yield  integral  representations  of  the  symmetric  and,  when  appropriate,  the  an¬ 
tisymmetric  part  of  D*.  Variational  formulations  of  the  effective  parameter  problem  for  D* 
are  given  in  [3,  15,  16]. 

Here  we  adapt  and  extend  the  mathematical  frameworks  developed  in  both  [3]  and  [41] 
to  the  case  of  a  time-independent  periodic  how.  We  also  discuss  how  the  frameworks  are 
extended  to  randomly  perturbed  periodic  hows.  In  particular,  for  each  approach  in  [3, 
41],  we  develop  Stieltjes  integral  representations  for  both  the  symmetric  and  antisymmetric 
parts  of  D*,  involving  the  molecular  diffusivity  and  a  spectral  measure  of  a  self-adjoint 
operator  that  depends  only  on  the  fluid  velocity  held.  Moreover,  for  each  approach  in  [3,  41], 
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we  provide  a  mathematical  foundation  for  rigorous  computation  of  D*  through  discrete, 
matrix  formulations  of  these  effective  parameter  problems.  For  an  alternative  computational 
method  that  is  accurate  for  large  Peclet  numbers  see  [20]. 

We  demonstrate  for  the  setting  introduced  in  [3],  that  the  discrete  spectral  measure 
is  given  explicitly  in  terms  of  the  eigenvalues  and  eigenvectors  of  a  standard  eigenvalue 
problem.  We  develop  a  projection  method  which  allows  the  spectral  measure  to  be  obtained 
from  the  eigenvalues  and  eigenvectors  of  a  much  smaller  matrix,  greatly  increasing  the 
efficiency  of  numerical  computations.  For  the  setting  described  in  [41],  we  show  that  the 
spectral  measure  is  instead  given  in  terms  of  the  eigenvalues  and  eigenvectors  of  a  generalized 
eigenvalue  problem. 

We  provide  a  detailed  matrix  analysis  that  demonstrates  the  two  approaches  yield  equiva¬ 
lent  spectral  representations  of  the  effective  diffusivity  tensor  D*  when  the  matrix  Laplacian, 
at  the  heart  of  both  approaches,  is  of  full  rank.  Moreover,  we  demonstrate  that  both  ap¬ 
proaches  can  be  formulated  in  terms  of  a  common,  standard  eigenvalue  problem.  However, 
the  matrix  Laplacian  becomes  rank  deficient  when  periodic  boundary  conditions  are  im¬ 
posed.  We  generalize  both  approaches  to  this  rank  deficient  setting  and  show  that  the  two 
generalized  approaches  again  yield  equivalent  spectral  representations  of  D*.  This  analysis 
provides  a  numerically  efficient  way  to  compute  the  discrete  spectral  measure  underlying  the 
integral  representation  for  D*.  We  utilize  this  unified  mathematical  framework  to  compute 
the  effective  diffusivity  for  some  model  flows. 


A.  Synopsis  of  the  paper 

In  Section  II,  we  formulate  the  effective  parameter  problem  for  advection  enhanced  diffu¬ 
sion  by  random,  incompressible  flows.  In  particular,  we  review  the  problem  of  homogenizing 
the  advection-diffusion  equation  for  such  flows,  developed  in  [31],  which  yields  a  rigorous 
functional  representation  for  the  effective  diffusivity  tensor  D*  in  terms  of  the  solution  of  a 
random  “cell  problem”  which  arises  in  the  homogenization  procedure. 

In  Section  III  A,  we  discuss  how  the  incompressibility  of  the  fluid  velocity  field  allows  it 
to  be  expressed  in  terms  of  the  divergence  of  an  antisymmetric  random  matrix.  This,  in 
turn,  allows  the  cell  problem  to  be  written  as  a  steady  state  diffusion  equation  involving  a 
curl-free  random  held.  Moreover,  we  demonstrate  how  this  leads  to  functional  formulas  for 
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the  symmetric  and  antisymmetric  parts  of  the  effective  diffusivity  tensor  D*,  involving  the 
curl-free  held  and  a  self-adjoint  random  operator.  In  Section  III  B,  we  discuss  how  the  cell 
problem  can  be  transformed  into  a  resolvent  formula  involving  this  self-adjoint  operator, 
which  acts  on  the  Hilbert  space  of  curl-free,  random  vector  fields  [2,  3].  We  then  discuss 
how  this  result  and  the  spectral  theorem  [46,  51]  yields  Sticltjes  integral  representations  for 
both  the  symmetric  [2,  3]  and  antisymmetric  parts  of  D*,  involving  a  spectral  measure  of 
the  operator.  We  utilize  analytical  properties  of  the  Stieltjes  integrals  to  provide  rigorous 
bounds  for  the  components  of  these  homogenized  tensors. 

In  Section  IV  A,  we  provide  a  mathematical  foundation  for  rigorous  computation  of  D* 
associated  with  this  approach.  In  particular,  a  discrete  representation  of  the  cell  problem 
involving  a  Hermitian  random  matrix  leads  to  Sticltjes  integral  representations  for  the  sym¬ 
metric  and  antisymmetric  parts  of  D*,  involving  a  discrete  spectral  measure  which  is  given 
explicitly  in  terms  of  the  eigenvalues  and  eigenvectors  of  the  matrix.  In  Section  IV  B,  we 
develop  a  projection  method  which  allows  spectral  statistics  to  be  obtained  by  diagonalizing 
a  much  smaller  random  matrix,  greatly  increasing  the  efficiency  of  numerical  computations. 

In  Section  V,  we  formulate  another  approach  to  the  effective  parameter  problem  [41], 
which  is  different  from  the  approach  discussed  in  Section  IV.  We  demonstrate  this  approach 
also  provides  Stieltjes  integral  representations  for  the  symmetric  and  antisymmetric  parts 
of  D*.  In  particular,  we  transform  the  cell  problem  into  a  resolvent  formula  involving  a 
self-adjoint  random  operator  acting  on  a  Sobolev  space  of  random  scalar  fields,  which  is  also 
a  Hilbert  space.  This  leads  to  functional  formulas  for  D*  and  a  resolvent  formula  for  the 
cell  problem.  This,  in  turn,  leads  to  Sticltjes  integral  representations  for  the  symmetric  and 
antisymmetric  parts  of  D*,  involving  a  spectral  measure  of  the  random  operator. 

The  symmetry  of  the  random  operator  depends  intimately  on  the  Sobolev-type  inner- 
product  in  this  approach.  Consequently,  its  matrix  formulation  is  substantially  different 
from  that  of  Section  IV.  This  technical  difficulty  is  overcome  in  Section  VI  by  casting  the 
effective  parameter  problem  in  terms  of  a  generalized  eigenvalue  problem,  which  has  the 
Sobolev-type  inner-product  as  a  key  feature.  This  leads  to  Stieltjes  integral  representations 
for  the  symmetric  and  antisymmetric  parts  of  D*,  involving  a  discrete  spectral  measure 
which  is  given  explicitly  in  terms  of  the  associated  generalized  eigenvalues  and  eigenvectors. 

The  inverse  Laplacian  operator  is  central  to  both  of  the  continuum  formulations  of  the 
effective  parameter  problems  described  in  Sections  III  and  V.  Consequently,  the  matrix 
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Laplacian  is  also  central  to  both  of  the  discrete  formulations  described  in  Sections  IV  and  VI, 
which  is  assumed  to  be  of  full  rank  so  that  it  is  invertible.  Given  this  condition,  we  demon¬ 
strate  in  Section  VII  that  the  two  discrete  formulations  of  the  effective  parameter  problems 
given  in  Sections  IV  and  VI  produce  equivalent  spectral  representations  for  the  effective 
diffusivity  tensor  D*. 

In  Section  VIII,  we  generalize  the  mathematical  frameworks  formulated  in  Sections  IV 
and  VI  to  the  case  of  a  rank  deficient  matrix  Laplacian.  This  is  important  in  computations 
of  the  effective  diffusivity  tensor  D*  for  randomly  perturbed  periodic  flows,  as  the  matrix 
Laplacian  with  periodic  boundary  conditions  is  rank  deficient.  A  detailed  matrix  analysis 
shows  that  both  approaches  yield  equivalent  spectral  representations  of  D*  and  that  the 
spectral  statistics  can  be  obtained  from  a  standard  eigenvalue  problem  of  a  much  smaller 
random  matrix.  When  the  matrix  Laplacian  is  of  full  rank,  these  generalized  formulations 
reduce  to  the  formulations  in  Sections  IV  and  VI. 

In  Section  IX  we  numerically  compute  the  spectral  measures  and  effective  diffusivity 
tensors  for  several  non-random  periodic  flows  as  well  as  randomly  perturbed  periodic  flows. 
Our  computations  are  in  excellent  agreement  with  theoretical  results  [3]  for  shear  flow.  We 
relate  spectral  characteristics  of  our  computations  to  flow  geometry  and  transport  properties. 

II.  HOMOGENIZATION  OF  THE  ADVECTION-DIFFUSION  EQUATION 

The  dispersion  of  a  cloud  of  passive  scalars  with  density  0  diffusing  with  molecular  dif¬ 
fusivity  £  and  being  advected  by  a  incompressible  velocity  held  u  satisfying  V-u  =  0  is 
described  by  the  advection-diffusion  equation 

4>t(t,  x)  =  u(x) •'V cj)(t ,  x)  +  eA(j)(t,  x),  0(0,  x)  =  <fro{x),  x  e  t  >  0,  (1) 

with  initial  density  0o(a?)  given.  Here,  04  denotes  partial  differentiation  of  0  with  respect 
to  time  t,  A  =  V-V  =  V2  is  the  Laplacian,  £  >  0,  d  is  the  system  dimension,  and  we 
denote  by  0  the  null  element  on  all  linear  spaces  in  question.  Moreover,  and  f 

is  the  operation  of  complex-conjugate-transpose,  with  =  |£|2.  Later,  we  will  extensively 
use  this  form  of  the  dot  product  over  complex  fields,  with  built  in  complex  conjugation. 
However,  we  stress  that  all  quantities  considered  in  this  section  are  real-valued.  We  assume 
for  now  that  the  time-independent  fluid  velocity  held  u  is  spatially  periodic  on  the  region 


V  C  Wl.  Later,  we  will  discuss  the  case  of  an  array  of  randomly  perturbed,  periodic  flows. 

The  long  time,  large  scale  dispersion  of  the  passive  scalars  can  be  described  [54]  by  an 
effective  diffusivity  tensor  D*.  An  explicit  representation  for  D*  can  be  found  using  methods 
of  homogenization  theory  [6,  39].  These  methods  have  demonstrated  that  the  averaged  or 
homogenized  behavior  of  the  advection-diffusion  equation  in  (1)  is  determined  by  a  diffusion 
equation  involving  an  averaged  scalar  density  0  and  a  (constant)  effective  diffusivity  tensor 
D*  [31] 


<k{t,x)  =  V-[D*V0(t,aj)],  <j>(0,x)  =  <j>o(x).  (2) 

The  components  D*k,  j,  k  =  1, . . . ,  d,  of  D*  are  given  by  [31] 

Djk  =  £Sik  +  (ujXk),  (3) 

where  (•)  denotes  volume  averaging  over  the  period  cell  V.  The  function  Xj  in  (3)  satisfies 
a  cell  problem  which  is  a  steady  state  advection-diffusion  equation  with  a  forcing  term 
involving  Uj ,  the  jtli  component  of  the  fluid  velocity  held  u  [16,  31], 

u-Vxj  +  eA  xj  =  - uj ,  (Vxj)  =  0.  (4) 

Equations  (2)-(4)  follow  from  the  assumption  that  the  initial  density  0o  is  slowly  varying 
relative  to  the  variations  of  the  velocity  held  u  [16,  31].  This  information  is  incorporated 
into  equation  (1)  by  introducing  a  small  dimensionless  parameter  5  <C  1  and  writing  [31] 

0(0,  x)  =  (j)0(Sx).  (5) 

Anticipating  that  0  will  have  diffusive  dynamics  as  t  — >  oo,  space  and  time  are  rescaled 
by  x  i — *  x/5  and  t  h- >  t/S2.  As  6  — >  0,  the  associated  solution  <ps(t,x )  =  0(t/h2,  x/5)  of 
equation  (1)  (in  the  rescaled  variables)  converges  to  0(i,  x)  which  satishes  equation  (2).  The 
convergence  is  in  an  L2  sense  that  depends  on  the  technical  assumptions  made  about  the 
fluid  velocity  held  u  [3,  15,  16,  29,  31,  41]. 

We  stress  that  the  cell  problem  in  (4)  involves  only  the  fast  variables  x/5  and  t/52  [31]. 
Other  space-time  scalings  have  also  been  considered,  which  have  lead  to  space-time  depen¬ 
dent  D*  [41]  and  even  anomalous  diffusive  dynamics  [29].  Homogenization  theorems  for 
space-time  dependent  fluid  velocity  helds  are  treated  in  [9,  29,  41]. 
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In  our  analysis  of  the  effective  diffusivity  tensor  D*,  it  will  be  convenient  to  use  non- 
dimensional  parameters.  We  therefore  assume  that  equation  (1)  has  been  non-dimensionalized 
as  follows.  Let  £  and  r  be  typical  length  and  time  scales  associated  with  the  problem  of 
interest.  Mapping  to  the  non-dimensional  variables  1 i— *•  t/r  and  x  i— ►  x/£,  one  finds  that  (j) 
satisfies  the  advection-diffusion  equation  in  (1)  with  a  non-dimensional  molecular  diffusivity 
and  fluid  velocity  held, 

£^Te/£2,  U  I — *  t  u/ £.  (6) 


In  the  case  of  a  time-independent,  spatially  periodic  how,  a  natural  choice  for  £  and  r  is, 
respectively,  the  maximum  cell  period  and  r  =  £/(\u\ 2)7/2,  so  that  (6)  is  given  by  i — ^ 
e/(i  ( | "ix | 2) 1//2)  and  u  >  u/ (\u\2)1/'2 .  In  this  case,  a  natural  definition  of  the  Peclet  number 
Pe  is 


Pe  = 


£{\u\2)1/2 


so  that  the  scaled  molecular  diffusivity  satisfies  £  =  Pe 


-l 


(7) 


III.  CURL-FREE  FIELDS  AND  THE  EFFECTIVE  DIFFUSIVITY  TENSOR 

In  this  section,  we  adapt  and  extend  a  method  [2,  3]  which  leads  to  integral  representations 
for  the  effective  diffusivity  tensor  D*.  More  specifically,  in  Section  III  A  we  provide  functional 
formulas  for  the  symmetric  and  antisymmetric  parts  of  D*,  involving  the  curl-free  vector 
held  V Xj  in  equation  (4) .  We  review  a  Hilbert  space  formulation  of  this  effective  parameter 
problem  [2,  3,  15,  16]  in  Section  III  B,  which  leads  to  a  resolvent  formula  for  Vyq,  involving 
a  self-adjoint  operator.  Moreover,  we  use  this  result  and  the  spectral  theorem  [46,  51]  to 
provide  Stieltjes  integral  representations  for  both  the  symmetric  and  antisymmetric  parts  of 
D*,  involving  a  spectral  measure  of  the  operator. 


A.  Functional  formulas  for  the  effective  diffusivity  tensor 


Since  u(x)  is  incompressible,  V*tt  =  0,  there  is  a  real  (non-dimensional)  antisymmetric 
matrix  H(:r)  [2,  3]  such  that 


H7  = 


-H, 


u  =  V-H, 


(8) 
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where  H7  denotes  transposition  of  the  matrix  H.  Due  to  the  anti-symmetry  of  the  matrix  H 
and  the  symmetry  of  the  Hessian  operator  VV  when  acting  on  a  sufficiently  smooth  space 
of  functions,  we  have  that  H  :  VVy?  =  0  for  all  such  smooth  functions  tp.  Consequently, 
V-[HVy?]  =  [V-H]-Vy9  +  H  :  VVy?  =  [V-H]-Vy>,  or  V-[HV]  =  [V-H]-V  as  operators 
acting  on  such  functions.  Using  this  identity  and  the  representation  of  the  fluid  velocity 
held  u  in  (8),  equation  (1)  can  be  written  as  a  diffusion  equation, 

0t  =  V-[DV0],  0(O,®)  =  0o(*),  D  =  el  +  H.  (9) 

Moreover,  the  cell  problem  in  (4)  can  be  written  as  a  steady-state  diffusion  equation  [15,  16], 

V-[D(Vyfc  +  efc)]  =  0,  (Vyfc)  =  0,  k  =  l,...,d.  (10) 

Here,  (•)  denotes  volume  averaging  over  the  period  cell  V,  e/,:  is  a  standard  basis  vector, 
k  =  1  and  D  (t,x)  =  e\  +  H(t,  a?)  can  be  viewed  as  a  local  diffusivity  tensor  with 

coefficients 


Djk  &&jk  Hjfc)  J)  k  1,  •  •  • ,  d,  (11) 

where  8jk  is  the  Kronecker  delta  and  I  is  the  identity  operator  on  Wl. 

The  symmetric  S*  and  antisymmetric  A*  parts  of  the  effective  diffusivity  tensor  D*  are 
defined  by 

D*  =  S*  +  A*,  S*  =  i  (D*  +  [D*]T)  ,  A*  =  i  (D*  -  [D*]T)  .  (12) 

Substituting  into  equation  (3)  the  expression  for  Uj  in  (4)  and  using  equation  (8),  u  =  V-hh 
one  finds  that  the  components  S*k  and  A*k,  j,k  =  1 , ,d,  of  S*  and  A*  can  be  written  in 
terms  of  the  following  functionals  involving  the  real-valued  vector  held  Vyy- 

s*„  =  Cfe  +  (VxrVxt)),  A-t  =  (H  VxrVxk).  (13) 

The  symmetry  S*fc  =  S *k-  of  the  tensor  S*  in  (13)  follows  from  the  fact  that  the  vector  held 
'Vxk  is  real-valued  so  that  (V^y- Vy^)  =  (Vxk’^Xj)-  Moreover,  the  positivity  condition 
(Vy'fc-Vyfc)  =  (|Vyfe|2)  >  0  demonstrates  that  the  effective  transport  of  the  scalar  density 
cf)  is  always  enhanced  in  all  of  the  principle  directions  e*,,  k  —  1, . . . ,  d,  by  the  presence  of 
an  incompressible  velocity  held,  Dkk  =  S*kk  >  e.  The  equality  Dkk  =  S*kk  follows  from  the 
skew-symmetry  of  the  matrix  A*,  A*kj  =  ~A*k,  hence  A kk  =  0.  This,  in  turn,  follows  from 
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the  skew-symmetry  of  the  real-valued  matrix  H,  which  implies  that  A*k  =  (HVx-pVxfc)  = 
— (Vxj-HV**)  =  -(HVxt-Vxy  =  -A’kj . 

We  conclude  this  section  by  noting  that  the  cell  problem  in  equation  (10)  is  equivalent  [2, 
3,  15,  16]  to  the  quasi-static  limit  of  Maxwell’s  equations  [19,  24,  35],  which  describe  the 
transport  properties  of  an  electromagnetic  wave  in  a  composite  material, 

V  X-Efc  =  0,  V-J7  =  0,  Jk  =  D-Efc,  (-Eh)  —  ek,  D  =  e\  +  H.  (14) 

Here,  Ek  =  Vy*,  +  e*,  plays  the  role  of  the  local  electric  held,  J k  =  D Ek  plays  the  role 
of  the  local  current  density,  and  D  =  el  +  H  plays  the  role  of  the  local  conductivity  tensor 
of  the  medium.  Since  H  is  skew-symmetric,  the  intensity-flux  relation  Jk  =  D Ek  is  not 
the  usual  Fourier  law,  but  instead  resembles  that  of  a  Hall  medium  [15,  16,  23,  35].  In 
Section  III  B,  we  employ  the  representation  of  the  cell  problem  in  (10)  and  adapt  the  analytic 
continuation  method  [19]  for  characterizing  transport  in  composites  to  provide  Stieltjes 
integral  representations  for  both  the  symmetric  and  antisymmetric  parts  of  the  effective 
diffusivity  tensor  D*,  involving  a  spectral  measure  of  a  self-adjoint  operator. 

B.  Hilbert  space  and  integral  representations 

The  analytic  continuation  method  provides  Stieltjes  integral  representations  for  the  bulk 
transport  coefficients  of  composite  media  [19] .  This  method  is  based  on  the  spectral  theorem 
of  Hilbert  space  theory  and  a  resolvent  formula  for,  say,  the  electric  held,  involving  a  self- 
adjoint  operator  [19]  or  matrix  [37]  which  depends  only  on  the  composite  geometry.  In  this 
section,  we  adapt  this  method  [2,  3]  to  provide  Stieltjes  integral  representations  for  both 
the  symmetric  and  antisymmetric  parts  of  the  effective  diffusivity  tensor  D*,  which  encodes 
the  complicated  geometry  of  the  fluid  velocity  held  in  a  spectral  measure  of  a  self-adjoint 
operator. 

The  analytical  platform  of  the  analytic  continuation  method  is  a  Hilbert  space  .  In 
the  effective  parameter  problem  for  effective  diffusivity,  the  mathematical  structure  of 
depends  on  the  specihc  details  of  the  fluid  velocity  held  of  interest.  When  one  considers  a 
fluid  velocity  held  that  is  spatially  periodic  on  a  region  V  C  Wl,  the  Hilbert  space  Jt?  can 
be  taken  to  be  [29] 


J4?  =  {$,  G  ®j[=1L2(V,  m)  :  £(x)  is  periodic  in  V}, 


(15) 
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where  m  is  the  normalized  Lebesgue  measure  (uniform  distribution)  on  restricted  to 
V,  and  L2(V,m)  is  the  space  of  complex- valued  scalar  functions  that  are  square-integrable 
with  respect  to  m.  The  Hilbert  space  is  equipped  with  a  sesquilinear  inner-product  (•,  •) 
defined  by  (£,£)  =  (£<),  with  «,£)  =  (£,  C)  for  £,  C  G  which  induces  a  norm  ||  ■  || 
defined  by  ||£||  =  (£,£)fo2  and  £,  G  Jf?  implies  that  ||£||  <  oo.  Here,  (•)  denotes  volume 
average  over  the  period  cell  V  with  respect  to  the  measure  m,  a  is  the  complex-conjugate 
of  the  scalar  a,  and  we  stress  that  the  dot  product  includes  the  operation  f  of 

complex-conjugate-transpose. 

One  could  also  imagine  a  random  fluid  flow  filling  all  of  Wl,  with  a  velocity  held  u  de¬ 
termined  by  the  probability  space  (0,  P )  with  cr-algebra  generated  by  the  sets  {u(x)  g5}, 
where  x  G  and  B  is  a  Borel  subset  of  Wl  [3].  Here,  fl  is  the  set  of  all  geometric  realiza¬ 
tions  of  u ,  which  is  indexed  by  the  parameter  w  G  h  representing  one  particular  geometric 
realization,  and  P  is  the  associated  probability  measure.  The  underlying  Hilbert  space  in 
this  case  can  be  taken  to  be  Jt?  =  L2( 0,  P ),  i.e.,  the  space  of  all  P-measurable  vector- valued 
functions  $,  satisfying  ||^||  =  (|C|2)1//2  <  oo,  where  (•)  denotes  ensemble  averaging  and  the 
underlying  sesquilinear  inner-product  is  defined  by  (£,  £)  =  (£•£)  .  In  this  case,  one  could 
consider  a  random  fluid  how  with  a  velocity  held  u  that  is  stationary  [31]  or  ergodic  [2,  3], 
for  example,  with  regularity  conditions  at  infinity,  i.e.,  as  |a?|  — >  oo.  In  these  cases,  one  works 
with  an  infinite  medium  directly,  which  presents  substantial  computational  difficulties. 

A  more  computationally  tractable,  random  system  is  given  by  a  n  x  n  array  of  randomly 
perturbed  periodic  hows  [16].  In  this  case,  the  cr-algebra  associated  with  the  underlying 
probability  space  is  generated  by  the  Lebesgue  measurable  subsets  of  Wl.  Here,  the  Hilbert 
space  is  given  by  equation  (15)  and  averaged  quantities  depend  on  the  realization  of  the 
random  medium  because  (•)  is  given  by  volume  averaging  over  the  period  cell  V  [16].  The 
effective  diffusivity  tensor  D*  is  obtained  by  taking  an  infinite  volume  limit,  D*  =  1  im,woc  D* , 
of  the  finite  volume  effective  diffusivity  tensor  D*  and  evoking  an  ergodic  theorem  [16,  19]. 
Numerically,  it  is  natural  to  spatially  average  each  statistical  trial  and  then  ensemble  average 
over  all  of  the  sampled  random  realizations. 

In  any  case,  once  the  Hilbert  space  Jt?  is  established,  with  associated  average  (•),  inner- 
product  (•,  •),  and  norm  ||  •  ||,  the  spectral  theory  presented  in  the  remainder  of  this  section 
progresses  independent  of  the  underlying  details,  as  it  lays  on  an  axiomatic  foundation  [51]. 
For  the  sake  of  numerical  tractability,  we  will  assume  that  the  Hilbert  space  J#7  is  given  by 
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equation  (15).  The  fluid  velocity  field  u  can  be  assumed  to  represent  a  periodic  or  randomly 
perturbed  periodic  flow. 

Now  consider  the  associated  Hilbert  space  J^x  of  curl-free  fields  [3,  15,  16,  19,  35], 

J*i?x  =  {£  G  :  Vx^  =  0  weakly  and  (£)  =  0}  .  (16) 

The  curl-free  vector  field  in  the  cell  problem  in  (10)  is  mean-zero  (Vy*,)  =  0.  When 
the  matrix  D  in  equation  (9)  is  bounded  in  the  operator  norm  ||  •  ||  induced  by  the  ^-inner- 
product  [18],  ||D||  <  oo,  then  there  exists  unique  G  J^fx  satisfying  equation  (10)  [19,  39]. 
We  assume  that 


0  <  £  <  oo,  ||H||  <  oo,  (IT) 

which  together  imply  that  ||D||  <  oo. 

The  linear  operator  T  =  V(A_1)V-  is  a  projection  onto  the  Hilbert  space  M’x  in  the 
sense  that  T  :  Jif  i— >•  fflx  and  (weakly)  for  all  $,  G  in  particular  T  Vy*,  =  Vy*,. 

It  is  based  on  convolution  with  respect  to  the  Green’s  function  for  the  Laplacian  A  [34,  49] . 
Applying  the  integro-differential  operator  VA-1  to  the  cell  problem  in  equation  (10)  yields 
r[(el  +  H)(Vyfc  +  ek)}  =  0.  Since  Tek  =  0  and  TVyfc  =  Vy fc,  this  is  equivalent  to 
(el  +  THr)  Vy k  =  —  THefc,  which  yields  the  following  resolvent  formula  for  Vy k 

VXk=(e\  +  A)~1gk,  A  =  THT,  gk  =  -THek.  (18) 

Since  T  is  a  projection  operator  onto  Jifx  C  Jf?  it  is  bounded  by  unity  in  operator  norm  on 
||r||  <  1  [18,  47].  Integration  by  parts  and  the  symmetry  of  the  operator  (—A)'1  [49] 
shows  that  T  is  also  a  symmetric  operator,  i.e.,  (r£*£)  =  (£•]?£)  for  all  £  G  dff .  These  two 
properties  of  T  show  that  it  is  a  self-adjoint  operator  on  J4?  [46].  Using  this  property  and 
rvyfc  =  Vy fc,  we  may  write  A*fc  in  equation  (13)  as  A*k  =  (A  VyyVyfc).  Consequently, 
substituting  the  resolvent  formula  for  Vy k  in  equation  (18)  into  the  functional  formulas  in 
equation  (13)  yields 

$*jk  —  £($jk  +  ((£l  +  A)  1gj-{e  I  +  A)  1gk)),  (19) 

^*jk  —  (A  (si  +  A)  lgj-(e\Jr  A)  1gk). 

Since  T  is  self-adjoint  on  Jif,  the  anti-symmetry  of  the  matrix  H  implies  that  A  =  THT 
is  an  antisymmetric  operator  on  i.e.,  (A .£•£)  =  —  (£-A£).  We  stress  that  the  operator 
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A  depends  only  on  the  fluid  velocity  field  via  equation  (8).  By  equation  (17),  the  operator 
A  is  bounded  on  with  ||A||  <  ||H||  <  oo.  This,  the  skew-symmetry  of  A,  and  the 
sesquilinearity  of  the  ^-inner-product  imply  that  M  =  —zA,  where  i  =  \/r—l,  is  a  bounded 
symmetric  operator,  hence  self-adjoint  on  Jt?  with  ||M||  =  ||  A||  <  oo.  The  spectrum  E  of  the 
self-adjoint  operator  M  is  real- valued  with  spectral  radius  equal  to  its  operator  norm  [46], 
he., 

SC[— ||H||,||H||].  (20) 

The  spectral  theorem  for  bounded  linear  self-adjoint  operators  in  Hilbert  space  states  that 
there  is  a  one-to-one  correspondence  between  the  operator  M  and  a  family  of  self-adjoint 
projection  operators  {Q(A)}a££  —  the  resolution  of  the  identity  —  satisfying  [51] 

lim  Q(A)  =  0,  lim  Q(A)  =  I.  (21) 

A— >  ini  Zj  A — >  sup  2j 

Furthermore,  for  all  ^,(  £  ,  the  complex-valued  function  of  the  spectral  variable  A 

defined  by  /^(A)  =  (Q(A)£,£)  has  real  and  imaginary  parts  that  are  strictly  increasing 
and  of  bounded  variation  [51].  Therefore,  there  is  a  complex  Radon-Stieltjes  measure 
associated  with  p^( A)  [18,  50,  51].  The  spectral  theorem  also  states  that,  for  all  complex¬ 
valued  functions  /  G  L2(p^)  and  h  G  L2(p^),  there  exist  linear  operators  denoted  by  /( M) 
and  h(M)  which  are  defined  in  terms  of  the  sesquilinear  functional  (/( M)  £,  h(M)  (,  )  [51]. 
In  particular,  this  functional  has  the  following  integral  representation  involving  the  Stieltjes 
measure  p^ ,  for  all  £,  £  G  J4?x, 

/OO 

7(A)  h( A)  d/i€C(A),  Mec(A)  =  (Q(A)£,  C),  (22) 

-oo 

where  the  integration  is  over  the  spectrum  E  of  M  [46,  51]  and  /  denotes  complex  conjugation 
of  the  scalar  function  /. 

Since  a  Stieltjes  measure  v  has  the  property  [51]  f'‘  dis(X)  =  v(b )  —  z/(a),  equation  (21) 
implies  that  the  mass  /j^-  of  the  measure  p^  is  given  by 

/OO  /*  OO 

d/i«(A)=/  d(Q(AK,C)  =  (^C),  (23) 

which  is  bounded  in  the  sense  that  j/i^J  <  ||£||  ||C||  <  00  f°r  all  £,  C,  G  Due  to 

the  sesquilinearity  of  the  inner-product  and  the  fact  that  the  projection  operator  Q(A) 
is  self-adjoint  on  the  complex- valued  function  p^(X)  satisfies  p-^(A)  =  /%( A)  and 
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=  ||Q(A)£||2.  Consider  the  associated  positive  Sticltjes  measure  /i^  and  the  real¬ 
valued  functions 

Re  hcc(A)  =  \  (hct(A)  +  7%  (A))  ,  Im  p,cc(A)  =  (/xK(A)  -  ;%(A))  ,  (24) 

with  associated  signed  Stieltjes  measures  Re  /j^  and  Im  /i^  [18]. 

For  notational  simplicity  we  denote  by  Pjk(X)  =  (Q(A )gj,gk)  instead  of  pgjgk(X).  We  now 
demonstrate  that  the  spectral  theorem  in  (22)  provides  Sticltjes  integral  representations  for 
the  functional  formulas  in  equation  (19),  involving  a  spectral  measure  fijk  associated  with  the 
function  g,jk( A),  where  the  real- valued  vector  held  gk  =  —  THe^.  is  defined  in  equation  (18). 
From  equation  (23),  the  mass  of  the  measure  pjk  is  real- valued  and  satisfies 

4  =  (9j,9k)  =  (RHe,  •  THe,,)  =  (H^He,  •  cfc>,  \p°jk\  <  ||H||2  <  oo,  (25) 


where  we  have  used  that  T  is  a  self-adjoint  projection  operator  on  J4?x.  In  equation  (22), 
set  M  =  —  tA,  £  =  g},  and  C  —  9k-  Moreover,  for  the  hrst  formula  in  equation  (19) 
set  /(A)  =  h( A)  =  (e  +  ^A)^1,  and  in  the  second  formula  set  /(A)  =  i\(e  +  «A)_1  and 
h( A)  =  (e  +  *A)_1,  with  A  e  K  and  |A|  <  ||H||  <  oo.  Since  (e2  +  A2)-1  <  1/e2  <  oo  and 
A2(e2  +  A2)-1  <  1  for  all  0  <  £  <  oo,  it  is  clear  from  equation  (25)  that  the  functions  /  and 
h  defined  above  satisfy  /,  h  G  L2(pkk)  for  all  k  —  1, ...  ,d  and  0  <  £  <  oo.  Consequently,  the 
spectral  theorem  in  (22)  implies  that  the  functional  formulas  for  S*fc  and  A*k  in  equation  (19) 
have  the  following  Stieltjes  integral  representations 


sh  = e 


djk  + 


r  dMA)A 

Loo  e2  +  \2) 


A*  - 
- 


Ld/ijfe(A) 

./-oo  £2  +  A2 


(26) 


which  involve  the  complex  measure  pjk- 

We  now  show  how  the  integrals  in  (26)  for  S*k  and  A*k  can  be  represented  in  terms  of  the 
signed  measures  Re  / ijk  and  Im  pjk,  respectively,  associated  with  the  functions  in  (24).  Since 
Vyfc  and  H  in  (13)  are  real-valued,  we  have  from  (18)  the  following  symmetry  conditions 


((el  +  A)  ^-(el  +  A)  1gk)  —  ((el  +  A)  1flffe*(el  +  A)  1gj)  (27) 

(A(el  +  A)  1g^{e\  +  A)  1gk)  =  ((el  +  A)  1gk'A  (el  +  A)  1gj). 


These  symmetries,  the  sesquilinearity  of  the  MJ -inner-product ,  the  linearity  [51]  of  the  Sticlt¬ 
jes  integral  in  (22)  with  respect  to  the  function  /i^(A)  and  equation  (26)  yield 


Sjk  ~  £  $ik  + 


>jk 


dRe  gjk( A)  A 
e2  +  A2  )  ’ 


A*  — 


A  dim  /ijk{ A) 
e2  +  A2 


(28) 
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where  we  have  used  Re  =  (/iJ-jk(A)+/IJ-fc(A))/2  and  Im  /^(A)  =  (/^(A)  —  Jijk(\))/(2i). 

The  formulas  for  S*fc  and  A*k  in  (28)  were  computed  with  respect  to  the  standard  basis 
{ej},  through  the  definition  of  /ijfc(A)  =  (Q(A)grJ-, gk)  with  gk  =  —  THe^.  We  now  show 
that,  given  S*k  and  A*k,  j,k  =  1  the  effective  diffusivity  tensor  can  be  computed 

relative  to  any  directions.  This  is  due  to  the  bilinearity  of  the  inner-product  underlying  the 
definition  of  pjk( A).  More  specifically,  if  £,  £  G  Wl  are  arbitrary  directions  of  interest,  then 

(Q(A)TH^.  THC)  =  Thjk  ajh(QWgj,9k),  where  the  constants  dj  and  bk,  j,  k  —  1, _ ,  d,  are 

the  coordinates  of  the  vectors  $,  and  (,  with  respect  to  the  standard  basis.  This  immediately 
leads  to  integral  representations  for  the  effective  diffusivity  tensor  relative  to  any  desired 
directions. 

We  conclude  this  section  with  a  discussion  regarding  some  rigorous  bounds  for  S*fc  and 
A*k,  j,k  =  1 , ,d,  that  follow  from  the  integral  representations  in  (28).  We  assume  that 
0  <  £  <  cx)  throughout  our  discussion.  From  equation  (20),  the  spectrum  E  of  the  operator 
M  =  —  iA  is  a  bounded  subset  of  M.  Denote  A+  =  supE  and  |A|+  =  supAeS  |A|,  and  recall 
that  inf  age  A2  =  0  as  A  =  0  is  a  limit  point  [49]  of  the  compact  operator  M  [8,  29].  Since  fikk 
is  a  positive  measure  with  finite  mass  gkk,  the  inequalities  l/(e2  +  \2+)  <  l/(e2  +  A2)  <  1/e2 , 
holding  for  all  A  G  E,  yield  [18] 

£[1  +  Akk/(£2  +  ^+)]  <  $*kk  <  £  [1  +  Akk/z2]-  (29) 

It  may  be  that  fikk  =  0,  hence  S *kk  =  e,  e.g.,  shear  flow  orthogonal  to  the  kth  direction  [3,  15]. 

When  j  jtz  k,  Re  fijk  is  a  signed  measure.  There  are  unique,  positive  measures  Re  /j^, 
and  Re  pjk  such  that  Re  /ijk  =  Re  p+k  —  Re  pjk  [18].  Moreover,  associated  with  the  signed 
measure  Re  /ijk  is  its  total  variation  |Re  g\ jk  [18] 

Re  Ajk  =  Re  fijk  -  Re  Ajki  lRe  a\ jk  =  Re  A%  +  Re  Ajk-  (30) 

From  equation  (25)  the  measures  Re  n+k  and  Re  pjk  have  bounded  mass,  [Re  Ap,]°  and 
[Re  respectively,  thus  the  mass  |Re  g\°-k  of  the  measure  |Re  /i\ jk  is  also  bounded.  Since 
|S*fc|  <  J  d|Re  /i\jk(\)/(e2  +  A2)  [18],  the  upper  bound  in  equation  (29)  with  gkk  replaced  by 
| Re  g\°-k  holds  for  the  positive  quantity  S*k \ .  Our  numerical  results  in  Section  IX  indicate 
that  /ifcfc  =  | Re  n\ jk,  j  ^  k,  for  2D  flows  that  are  symmetric  about  the  line  y  —  x,  such  as  as 
cat’s  eye  flow,  yielding  the  bound  \S*k\  <  S*kk  which  holds  for  such  flows  (see  Figures  3  and  4 
below).  These  bounds  for  S*fc  can  be  improved  upon  by  separately  considering  the  positive 
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and  negative  contributions  of  the  integral  representation  for  S*k,  yielding 

[Re  /4fc]°  tRe  I** fe]°  /  c*  /  tRe  Ri-]°  tRe  %']° 

^  £2  +  X2+  £  ~  jk  ~  £  £  £2  +  XI  ’ 


(31) 


In  a  similar  way,  we  obtain  the  following  bounds  for  A*fc 


lAl  +  lImAl jk  ,  A*  .  |A|+  |Im/r|°fc  . 

—  S  S  - — - ,  J  T  K, 


£‘ 


£‘ 


(32) 


where  |Im/x|°fc  is  the  finite  mass  of  the  total  variation  |Im  fi\ jk  =  Im  /x^,  +  Im  / ijk  of  the  signed 
measure  Im  [ijk  =  Im  /i+k  —  Im  /i~k  [18] .  Bounds  on  can  also  be  obtained  using  variational 
methods  [3,  15,  16]  as  well  as  Pade  approximants  [3,  5]  of  Sticltjes  functions. 


IV.  DISCRETE  SETTING:  HILBERT  SPACE  OF  CURL-FREE  FIELDS 

Here,  we  provide  a  discrete  formulation  of  the  effective  parameter  problem  in  Section  III  B, 
which  will  be  used  in  Section  IX  to  compute  a  discrete  spectral  measure  jijk  involved  in 
discrete  versions  of  the  Sticltjes  integral  representations  for  the  effective  diffusivity  tensor 
D*  in  equation  (28).  More  specifically,  we  use  the  cell  problem  in  equations  (10)  and  (18), 
written  as  (H  +  A)Vyfc  =  flh,  to  express  the  discrete  spectral  measure  explicitly  in  terms 
of  the  eigenvalues  and  eigenvectors  of  a  matrix  representation  A  of  the  operator  A  =  THr, 
not  to  be  confused  with  the  antisymmetric  part  A*  of  the  effective  diffusivity  tensor.  In 
Section  IV  A,  we  discuss  the  spectral  properties  of  the  matrix  A  and  provide  the  Stieltjes 
integral  representations  for  the  effective  diffusivity  tensor  D*  in  (28),  involving  the  discrete 
spectral  measure  / ijk .  A  projection  method  is  formulated  in  Section  IV  B,  which  shows  that 
the  discrete  measure  fijk  is  determined  by  the  eigenvalues  and  eigenvectors  of  a  matrix  that 
is  much  smaller  than  the  matrix  A.  This  greatly  increases  the  efficiency  of  our  numerical 
computations  of  D*  in  Section  IX. 

A.  Matrix  formulation  of  the  effective  parameter  problem 

In  the  discrete  setting  [37],  H  is  represented  by  a  banded  antisymmetric  matrix.  For 
simplicity  we  will  not  make  a  notational  distinction  between  the  two  cases,  as  the  context 
will  be  made  clear.  The  differential  operator  V  is  represented  by  a  finite  difference  matrix 
V  [13,  37],  where  V7  =  (Vf, . . . ,  Vj)  and  Vj,  j  —  1, . . . ,  d,  are  also  finite  difference  matrices. 
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Moreover,  the  divergence  operator  V-  is  given  by  —  V7  and  the  matrix  representation  of 
the  negative  Laplacian  —A  is  given  by  V7  V.  Consequently,  the  projection  operator  T  = 
V(A_1)V-  is  represented  by  the  symmetric  projection  matrix  T  =  V(V7  V)_1VT,  satisfying 
T2  =  T  and  TV  =  V,  where  (V7  V)-1  is  now  interpreted  as  a  matrix  inversion.  We  assume 
here  that  the  matrix  V  is  of  full  rank  so  that  (V7V)_1  exists.  The  rank  deficient  case, 
where  the  matrix  V7  V  is  singular,  is  examined  in  detail  in  Section  VIII.  In  this  way, 
the  integro-differential  operator  A  is  represented  by  an  antisymmetric  matrix  A  =  THT 
satisfying  A7  =  —  A,  which  is  not  to  be  confused  with  the  antisymmetric  part  A*  of  the 
effective  diffusivity  tensor  D*.  In  a  similar  way,  the  vectors  gk  =  —  THe^,  k  —  1, . . . ,  d,  are 
redefined  in  this  matrix  setting  and,  for  simplicity,  we  will  not  make  a  notational  distinction 
between  the  two  cases  for  the  vectors  gk  and  ek,  as  the  context  will  be  made  clear. 

The  spectrum  E  of  the  antisymmetric  matrix  A,  of  size  N,  say,  consists  solely  of  eigen¬ 
values  vn,  n  =  1, . . . ,  N,  with  corresponding  eigenvectors  wn  satisfying  A wn  =  vnwn.  It 
is  well  known  [22]  that  its  eigenvalues  vn  are  purely  imaginary,  vn  =  i\n  with  An  G  M. 
Therefore,  the  matrix  M  =  —  zA  is  Hermitian  (IVb  =  M)  and  it  has  the  same  eigenvectors  wn 
as  the  matrix  A  and  real  eigenvalues  given  by  An  =  Im  vn .  It  is  also  well  known  [25]  that 
the  eigenvectors  wn,  n  —  1, . . .,  N,  form  an  orthonormal  basis  for  CN ,  i.e.,  w^wm  =  Snm 
and  for  every  £  G  CN  we  have  ^  =  J2n(wn^)wn  =  (Yln  wnwn)  £•  Consequently,  dehning 
Qn  =  wnw jj,  n  =  1, . . . ,  N,  to  be  the  mutually  orthogonal  projection  matrices  onto  the 
eigenspaces  spanned  by  the  wn , 

N 

^  ^  Q n  =  I;  Q n  =  QzQm  =  Q l  $lm-  (33) 

n=  1 

We  now  use  equation  (33)  to  prove  the  spectral  theorem  in  (22)  for  this  matrix  setting. 
From  =  A nwn  and  equation  (33)  we  have  that  MQn  =  AnQn.  This  formula  and  (33) 
then  imply  that  the  matrix  M  has  the  spectral  decomposition  M  =  Yn  YiQn-  By  the  mutual 
orthogonality  of  the  projection  matrices  Qn  and  by  induction,  we  have  that  Mm  =  Yln 
for  all  m  G  N.  This,  in  turn,  implies  that  /( M)  =  ]T)n  /( An)Qn  for  any  polynomial  /  :  M  i— *  C. 
From  the  mutual  orthogonality  and  the  symmetry  of  the  projection  matrices  Qn  it  follows 
that,  for  all  £  G  CN  and  complex- valued  polynomials  /(A)  and  h( A),  the  bilinear  functional 

(/( M)£»/i(M)£)  has  the  integral  representation  in  equation  (22),  with  M  substituted  by  M. 
Moreover,  the  complex-valued  function  =  (Q(A)^,C)  hi  equation  (22)  is  now  given 

by  p.^(A)  =  (Q(A)^-C),  where  the  associated  matrix  representation  Q(A)  of  the  projection 
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operator  Q(A)  and  the  discrete  spectral  measure  d/i^ ( A)  are  given  by 

Q(A)=  £  d(A  —  An)Qn,  d/x^(A)  —  ^  (£\„(dA)[Qn£‘£])-  (34) 

n:  An<A  n:  An<A 

Here,  0(A)  is  the  Heaviside  function,  satisfying  9( A)  =  0  for  A  <  0  and  0(A)  =  1  for  A  >  0, 
and  <5A„(dA)  is  the  5-measure  centered  at  An. 

The  spectral  theorem  for  symmetric  matrices  also  holds  for  functions  of  the  form  /(A)  = 
a\l(b  +  cA)~m,  with  /,  m  G  N,  a,  b,  c  G  C,  and  b  +  cA  ^  0  for  all  A  G  E.  We  will  demonstrate 
this  for  the  special  cases  that  arise  in  the  functional  formulas  for  S*k  and  A*fc  in  equation  (19), 
with  A  substituted  with  A,  which  yield  the  integral  representations  in  (26)  involving  a  dis¬ 
crete  spectral  measure  d/ij/c(A)  associated  with  the  measure  in  (34).  The  argument  involving 
the  function  /(A)  above  is  a  simple  extension  of  that  given  here. 

As  A  is  a  real- valued  matrix,  its  eigenvalues  vn  =  x\n,  n  =  1 ,N,  and  eigenvectors 
wn  come  in  complex-conjugate  pairs  [22],  Therefore,  if  the  size  A  of  A  is  even,  then  we 
may  re-number  the  index  set  In  as  In  =  { — A//2, . . . ,  — 1, 1, . . . ,  A/2}  such  that  V-n  =  nH  = 
—vn  and  iu_n  =  T<v  If  A  is  odd  then  vq  =  0  is  also  an  eigenvalue  with  a  real-valued 
eigenvector  Wq.  Denoting  by  W  the  matrix  with  columns  consisting  of  the  eigenvectors  wn, 
T  =  diag(i;_7v/2,  •  •  •  ,vn/ 2)  the  diagonal  matrix  with  eigenvalues  vn  on  the  main  diagonal, 
and  A  =  diag(A_jv/25  ■  •  • ,  A v/2),  we  have  that  A  =  WYV\A  =  i\ A/AW/  where  the  matrix  W 
is  unitary  V\AW  =  WV\A  =  I  [22],  Therefore,  the  matrix  M  =  —  tA  =  WAV\A  is  Hermit  ian 
(Mt  =  M). 

This  spectral  decomposition  of  A  demonstrates  that  the  matrix  (el  +  A)-1,  is  well  defined 
for  all  0  <  e  <  00.  In  particular,  since  V\A  =  W-1,  it  has  the  following  useful  representation 
(el  +  A)-1  =  W(el  +  zA)_1Wi,  where  (el  +  w\)-1  is  a  diagonal  matrix  with  entries  l/(e  +  *A). 
This  allows  the  discrete  version  of  the  resolvent  formula  in  equation  (18)  to  be  written  as 

VXj  =  W(el  +  *A)_1WtflrJ-,  9j  =  -THej.  (35) 

Substituting  the  resolvent  formula  for  VXk  in  equation  (35)  into  the  discrete  version  of  (13) 
and  using  TV  =  V  to  write  (H  VXj'^Xk)  =  (A  Vx^-Vx^,),  we  obtain  the  following  analogue 
of  equation  (19), 


SJt  =  e( Sjt  +  ((el  +  tA)  “Wtty.fcl+tA)  'W Vi)). 
A*t  =  (*A(eI  +  iA)-1W<g/-(el  + 


(36) 
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where  we  have  used  that  W'  =  W  1 .  The  quadratic  form  •  W'gA.  arising  in  (36)  can  be 
written  in  terms  of  the  projection  matrices  Qn  dehned  in  (33)  as  follows 

w  gk  =  (wigj){wigk)  =  Qngj-gk-  (37) 

ne/jv  n£lN 

This  implies  that  the  functional  formulas  for  S*k  and  A*fc  in  equation  (36)  have  the  in¬ 
tegral  representations  in  equation  (26)  with  discrete  spectral  measure  d/q-fc(A)  dehned  in 
equation  (34)  with  ^  =  g  •  =  — THe.,  and  C  =  9k  =  —  THe*,. 

As  in  the  abstract  Hilbert  space  setting  discussed  in  Section  III  B,  we  may  use  the  fact 
that  the  matrix  A,  the  vector  g  ■,  and  the  molecular  diffusivity  £  are  real- valued  to  obtain 
the  symmetry  relations  in  equation  (27),  with  A  changed  to  A.  This  allows  us  to  rewrite 
the  integral  representations  for  S*k  and  A*k  in  (26)  involving  the  discrete,  complex  measure 
dfj,jk(X),  as  the  integral  representations  in  equation  (28)  involving  the  real  signed  measures 
Re  jijk  and  Im  ji]k.  As  in  the  abstract  Hilbert  space  setting,  these  signed  measures  are 
determined  by  the  functions  Re  gjk{ A)  and  Im  /ijk( A)  in  equation  (24),  where  in  this  matrix 
setting  /j,jk( A)  =  (Q(A )g3-gk)  and  Q(A)  is  dehned  in  equation  (34).  As  the  projection  matrix 
Qn  is  Hermitian,  we  have  [Q ngk'9j\  =  [Q ngj-gk\.  Consequently,  from  equations  (24)  and  (34) 
we  have  that 

Re  Hjk{ A)  =  +  Qn)gj-gk])  (38) 

n:  An<A 

Im  Aifc(A)  =  7^-  (^(A-An)[(Q n-  Qn)gj-gk}), 

n:  \n<\ 

where  [(Qn  +  Q n)gfgk}  =  2Re  [Q n9j-gk]  and  [(Qn  -  Q n)gygk]  =  2/hn  [Q ngj-gk]- 

Since  the  eigenvalues  vn  =  i\n  and  eigenvectors  wn  of  the  matrix  A  come  in  complex 
conjugate  pairs,  the  representations  of  the  measures  Re  njk  and  Im  pjk,  which  follow  from 
the  functions  in  equation  (38),  can  be  simplified  and  shown  [41]  to  depend  only  on  the 
restricted  index  set  {n  >  0  :  \n  <  A}.  This  is  clear  from  equations  (28)  and  (38),  since  for 
n  >  0  we  have  X2_n  =  (— An)2  =  A2  and  w_n  =  wn,  thus  Q_n  =  Qn.  Consequently,  we  have 
that 

Re  [Q ngj-gk]  +  Re  \Q-ngj-gk\  =  2Re  \Qng3-gk]  (39) 

A„Im  [Q n9j'9k\  +  A-nlm  =  2AnIm  [Q„fly#fc], 


with  AoIm[Q0^-*^fc]  =  0. 
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B.  Projection  method 


We  now  formulate  a  projection  method  which  demonstrates  that  the  discrete  spectral 
measure  /ijk  in  (38)  is  determined  by  the  eigenvalues  and  eigenvectors  of  a  matrix  that  is 
much  smaller  than  the  matrix  A.  This  will  be  used  in  Section  IX  to  greatly  increases  the 
efficiency  of  our  numerical  computations  of  D*.  This  method  follows  from  the  projective 
nature  of  the  matrix  T,  which  causes  A  =  THT  to  have  a  large  null  space  associated  with 
the  null  space  of  T.  We  demonstrate  that  the  spectral  weights  [Qngj-gk]  associated  with 
this  null  space  are  identically  zero  and  therefore  do  not  need  to  be  computed  at  all. 

Since  T  is  a  real-symmetric  projection  matrix  of  size  N,  its  eigenvalues  yn,  n  —  1, . . . ,  N, 
satisfy  =  0, 1  [22],  Consequently,  T  has  the  spectral  decomposition  T  =  PGP7 ,  where  P  is 
an  orthogonal  matrix,  PP7  =  P7  P  =  I,  with  columns  consisting  of  the  eigenvectors  of  T,  and 
the  diagonal  matrix  G  has  the  eigenvalues  along  its  main  diagonal.  Write  P  =  [Po  Pi], 
where  the  columns  of  the  N  x  No  matrix  Po  and  the  N  x  N\  matrix  Pi  are  orthonormal 
eigenvectors  that  span  the  null  space  and  range  of  T,  respectively,  with  No  +  N\  =  N  [22],  ft 
follows  that  G  =  diag(Ov0,  IjvJj  where  Ojv0  and  1^  are  vectors  of  zeros  and  ones  of  length 
No  and  N\ ,  respectively.  This,  in  turn,  implies  that  the  matrix  T  can  be  written  as 

r  =  PiPf.  (40) 


The  matrix  P  being  orthogonal  implies  that  P7  Pi  =  In,  where  In  is  the  identity  matrix  of 
size  N\  x  N\.  This  demonstrates  that  the  matrix  T  satisfies  T2  =  T.  Moreover,  T  projects 
vectors  in  MiV  on  to  the  subspace  spanned  by  the  columns  of  Pi. 

Using  the  spectral  decomposition  T  =  PGP7  of  T,  we  may  write  the  matrix  A  =  THT  as 
A  =  P[G(P7  HP)G]P7  .  The  block  matrix  form  of  the  matrix  G(PrHP)G  =  P7  AP  is  given  by 


G(PtHP)G  = 


Ooo  O0i 
Oio  In 


Po  Pi 


Ooo  O0i 

Oio  In 


Ooo  O0i 

Oio  PfHPi 


where  0 ab  is  a  matrix  of  zeros  of  size  Na  x  Nb,  a,b  =  0,1. 

Due  to  the  skew-symmetry  of  H,  the  matrix  PfHPi  of  size  N\  x  N\  is  also  antisymmetric 
and  consequently  has  the  spectral  decomposition 


PfHP^tRuAnRj,, 


(42) 
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where  Rn  is  a  unitary  matrix,  R^Rn  =  RnR^  =  In,  and  An  is  a  real- valued  diagonal 
matrix.  Equations  (41)  and  (42)  allow  the  matrix  A  to  be  written  as 

A  =  A/VAWt,  W  =  PR,  R  =  lo°  °01  ,  A  =  °00  °01  ,  (43) 

Oio  Rn  Oio  An 

where  loo  is  the  identity  matrix  of  size  No  x  No.  Since  Rn  is  a  unitary  matrix,  the  matrix  R 
also  is.  Consequently,  since  P  is  orthogonal,  W  is  a  unitary  matrix  V\AW  =  WV\C  =  I. 

Equation  (43)  demonstrates  that  the  eigenvalues  i\n  of  the  matrix  A  are  zero  for  all 
n  =  1, . . . ,  jV0.  We  now  show  that  the  associated  spectral  weights  [Q ngj-gk]  are  also  zero 
for  all  n  =  Here,  g •  =  — THe^,  Qn  =  wnw^,  and  wn,  n  =  1,  are  the 

eigenvectors  of  A  which  comprise  the  columns  of  the  matrix  W.  From  equation  (43),  we 
see  that  W*  =  RfPT.  Since  T  =  PGPT  and  P  is  an  orthogonal  matrix,  this  implies  that 
WW  =  R^GP7  .  It  follows  from  the  block  forms  of  the  matrices  G,  P,  and  R  in  equations  (41) 


and  (43) ,  that  WW  is  given  by 


Wfr  =  RfGPT  = 


loo  Ooi 
Oio  Rn 


Ooo  Ooi 
Oio  bi 


RliPf 


where  Oqn  is  a  matrix  of  zeros  of  size  No  x  N.  It  follows  from  gJ  =  — THe.,  and  equation  (44) 
that  wJlgJ  =  0  for  all  n  —  1, . . .  ,N0.  This  and  equation  (37)  imply  that  [Qngj-gk]  =  0  for 


all  n  =  1, . . . ,  iVo,  as  claimed. 


In  Section  IX  we  will  employ  this  projection  method  to  compute  the  discrete  spectral  mea¬ 
sure  d/ijfc(A)  by  directly  computing  the  eigenvectors  and  eigenvalues  of  the  matrix  P7  HPn 
Our  analysis  demonstrates  that  the  discrete  spectral  measure  dfijk(X)  associated  with  equa¬ 
tion  (34)  does  not  depend  on  components  [Q n  =  1,...,N0,  and  therefore  do  not 
need  to  be  computed  at  all.  Moreover,  in  the  case  of  a  randomly  perturbed,  periodic  velocity 
held,  since  the  matrix  T  is  non-random,  the  spectral  decomposition  T  =  PGP7  needs  to  be 
computed  only  once,  while  the  spectral  decomposition  in  (42)  of  the  much  smaller  matrix 
P7  HPi  of  size  N\  x  N\  is  performed  repeatedly  to  gather  spectral  statistics.  This  greatly 
increases  the  efficiency  of  our  associated  computations. 


V.  SOBOLEV  SPACE  AND  THE  EFFECTIVE  DIFFUSIVITY  TENSOR 

In  this  section,  we  adapt  and  extend  an  alternate  method  [8,  41]  to  the  method  presented 
in  Section  III,  which  provides  integral  representations  for  the  effective  diffusivity  tensor 
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D*,  and  leads  to  a  more  direct  approach  for  its  computation  than  the  projection  method 
discussed  in  Section  IV  B.  We  provide  functional  formulas  that  are  analogous  to  the  formulas 
in  equation  (13)  for  the  symmetric  S*  and  antisymmetric  A*  parts  of  D*,  involving  the 
scalar  field  Xj  in  equation  (4) .  We  also  provide  a  Sobolev  space  formulation  of  the  effective 
parameter  problem  [41]  which  yields  a  resolvent  formula  for  \j  that  is  analogous  to  the 
formula  in  (18),  involving  a  self-adjoint  operator  that  depends  only  on  the  fluid  velocity 
field.  This  leads  to  the  Stieltjes  integral  representations  in  (28),  involving  a  spectral  measure 
of  the  operator. 

In  contrast  to  Section  III  B,  which  defined  the  effective  tensor  D*  in  terms  of  the  vector 
field  Vyy,  here  we  define  D*  in  terms  of  the  scalar  field  Xj-  Consider  the  Hilbert  space  H, 

H  =  {/  G  L2(V,  m)  :  f(x)  is  periodic  in  V},  (45) 

where  m  is  the  Lebesgue  measure  on  Wl,  restricted  to  V,  and  the  cr-algebra  associated  with 
the  underlying  probability  space  is  generated  by  the  Lebesgue  measurable  subsets  of  M,d.  The 
Hilbert  space  H  is  equipped  with  a  sesquilinear  inner  product  (•,  •)  defined  by  (/,  h)  —  ( f  h ), 
with  ( h ,  /)  =  (/,  h)  for  /,  h  G  H,  which  induces  a  norm  ||  ■  ||  defined  by  ||/||  =  (/,  /)1//2  and 
/  G  H  implies  ||/||  <  oo.  Now  consider  the  associated  Sobolev  space  7 T1,2,  which  itself  is  a 
Hilbert  space  [8], 

W1'2  =  {/£«:  II  /II  i,2  <  OO,  {/)  =  0},  ll/ll!,  2  =  (IV/iyy  (46) 

where  ||  •  Hi^  is  the  norm  induced  by  the  underlying  sesquilinear  inner-product  (•,  -)\y2  defined 
by  (f,h)  1,2  =  (V f-'Vh),  with  (h,f) i,2  =  (f,h) lj2. 

Recall  the  definition  of  the  components  D*k  =  eSjk  +  (ujXk),  j,  k  —  1, . . . ,  d,  of  the  effective 
diffusivity  tensor  D*  in  (3).  Rewrite  the  functional  (■ UjXk )  as  [41], 

(ujXk)  =  ([AA  ^ufixk)  =  -(VA  _1«j-Vxfc)  =  —  (A_1rij,  Xk)i,2,  (47) 

where  we  have  used  the  periodicity  of  the  functions  Uj  and  Xk ■  Substituting  into  equa¬ 
tion  (47)  the  expression  —  Uj  =  u-'Vxj  +  ^A Xj  for  —  Uj  in  (4)  yields  the  following  functional 
formulas  for  the  components  S*fc  and  A*fc,  j,k  =  1 , ,d,  of  the  symmetric  S*  and  antisym¬ 
metric  A*  parts  of  D*  defined  in  equation  (12), 


S*jk  =  £(Sik+  (Xj,Xk)  1,2), 


=  (■ Axj,Xk)i,2 , 


A  =  A_1[t4.v;. 


(48) 
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The  functional  formulas  in  (48)  are  analogous  to  the  functional  formulas  in  equation  (13). 
Dne  to  the  incompressibility  of  the  fluid  velocity  field,  V-u  =  0,  the  operator  A  is  anti¬ 
symmetric  on  741,2,  i.e.,  (Af,h)  1,2  =  —  (/,  Ah)i^  for  all  f,h  G  74 12  (see  equation  (Al)). 
We  stress  that  the  operator  A  is  independent  of  e.  Since  the  scalar  fields  Xj  and  Ayy  are 
real-valued ,  we  have  that  (xj,Xk)  1,2  =  (Xk,Xj)  1,2  and  A*kj  =  (Axk,Xj)  1,2  =  = 

—  (Ax j,Xk)  1,2  —  —  A*fc,  confirming  that  S*  is  symmetric  and  A*  is  antisymmetric. 

Applying  the  operator  A-1  to  both  sides  of  the  the  cell  problem  eA Xj  +  u-W\j  =  — Uj 
in  equation  (4)  yields  the  following  resolvent  formula  for  \j  involving  the  antisymmetric 
operator  A,  which  is  analogous  to  equation  (18), 

Xj  =(£  +  A)-1gj ,  9j  =  -A  ~xUj.  (49) 

Substituting  the  resolvent  formula  for  Xj  in  (49)  into  the  functionals  in  equation  (48)  yields 

Sjfc  =  £(5jk  +  ((£  +  A)  lgj ,  (e  +  A)  1gk)  1,2),  (50) 

A =  (A  (e  +  A)  ^ Qj ,  (e  +  A)  1(7fc)i,2, 

which  is  a  direct  analogue  of  equation  (19). 

Since  V  is  a  bounded  domain,  the  linear  operator  A_1  is  bounded  on  74  [49].  When  \u\ 1 
is  uniformly  bounded  on  the  period  cell  V, 

sup  \u(x)\2  <  00,  (51) 

xGV 

the  linear  operator  A  is  bounded  on  741'2,  with  (see  equations  (A2)  and  (A3)) 

||A||i,2<  [  II A^1!!  sup \u(x)\'2  ]1/2  <  (X).  (52) 

x£V 

All  of  the  fluid  velocity  fields  that  we  consider  in  Section  IX  satisfy  equation  (51).  More 
generally,  for  uk  G  74,  k  —  1, . . . ,  d,  the  operator  A  is  compact  on  741'2  [8],  hence  bounded  [49]. 
It  follows  that  M  =  —  1A  is  a  bounded  linear  operator  on  741,2  with  ||M||ii2  =  || A|| ii2  <  00, 
where  %  =  \J — 1.  Moreover,  the  skew-symmetry  of  A  and  the  sesquilinearity  of  the  7 41,2- 
inner-product  imply  that  M  is  also  symmetric,  (Mf,h) ij2  =  (f,Mh) ii2,  hence  self-adjoint 
on  741,2  [46].  The  spectrum  X  of  the  self-adjoint  operator  M  is  real-valued,  with  spectral 
radius  equal  to  its  operator  norm  [46],  i.e., 


SC  [-\\A\\lt 2,||A||1j2]. 


(53) 
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Exactly  as  in  Section  III  B,  direct  analogues  of  equations  (22)-(28)  lead  to  the  Stieltjes 
integral  representations  in  (28),  involving  a  Stieltjes  measure  pjk  associated  with  the  function 
of  the  spectral  variable  A  defined  by  njk( A)  =  (Q(X)gj,  gk)i,2,  where  gj  =  —A ~xUj  is  defined 
in  (49)  and  {<5(A)}ags  is  the  family  of  self-adjoint  projection  operators  that  is  in  one-to- 
one  correspondence  with  the  bounded  linear  self-adjoint  operator  M  on  the  Hilbert  space 
H1'2.  From  equation  (23),  the  mass  of  the  measure  fijk  is  given  by  n°-k  =  (gj,gk)  1,2  = 
(V A~1Uj-'V A~2Uk)  so  that 

rfk  =  <[(-A)_1«j]“fc>,  |/4l  <  IIA_1m/II  IKII  <  00.  (54) 

VI.  DISCRETE  SETTING:  SOBOLEV  SPACE  OF  SCALAR  FIELDS 

In  Section  IX,  we  consider  a  discrete  approximation  of  the  cell  problem  in  (4)  written 
as  (e  +  iM)xj  =  gj,  where  M  =  —  %A,  A  =  A_1[tt*V],  and  gj  =  —A ~2Uj,  as  defined  in 
equations  (48)  and  (49).  From  the  formula  u  —  V-H  in  equation  (8)  and  our  discussion 
in  Section  IV  A,  we  may  write  the  discrete,  matrix  representation  M  of  the  self-adjoint 
operator  M  =  A-1[— ?V-HV]  by  M  =  (— VrV)-1[— ?V7  HV].  As  in  Section  IV  A,  for  sim¬ 
plicity,  we  do  not  make  a  notational  distinction  for  the  matrix  H,  between  the  continuum 
and  discrete  settings  as  the  context  will  be  clear.  This  composition  of  the  Hermitian  matrix 
— zVTHV  and  the  real-symmetric  matrix  (— V7  V)'1  is  not  symmetric.  From  equation  (Al) 
we  see  that  the  symmetry  properties  of  the  operator  M  depend  intimately  on  the  inner- 
product  (/,  h) i,2  =  (V/*V/i)  of  the  underlying  Sobolev  space  7i1:2  defined  in  equation  (46). 
Therefore,  the  properties  of  this  inner-product  must  be  included  in  the  discrete  formulation 
of  the  effective  diffusivity  tensor  D*. 

In  this  section,  we  provide  a  discrete,  matrix  formulation  of  the  effective  parameter  prob¬ 
lem  introduced  in  Section  V,  which  involves  a  generalized  eigenvalue  problem  that  has 
the  Sobolev-type  inner-product  as  a  central  feature.  In  particular,  this  formulation  re¬ 
tains  the  key  properties  of  the  weak  form  of  the  eigenvalue  problem  (M(pn,  </?n)i,2  =  An. 
Namely,  the  operator  M  is  self-adjoint  with  respect  to  the  inner-product  (•,  defined  by 
{f,h)  1,2  =  (V/*V/i),  its  eigenfunctions  <pn  G  7 i1,2  are  orthonormal  {(pn,  <Pm)  1,2  =  5nm, 
n,m  =  1,2,3,...,  with  respect  to  the  inner-product  (-, -)i,2,  and  the  spectrum  S  of  M  is 
real  valued,  Eel.  Towards  this  goal,  consider  the  eigenvalue  problem  Mipn  =  \n(pn  in  the 
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following  form, 

-lVMV<pn  =  An  A  (fin-  (55) 

Equation  (55)  is  well  defined  for  Hjfc  G  C1(V)  and  <pn  G  C2(V),  where  Cr(V)  the  space  of 
continuously  differentiable  functions  of  order  r  with  domain  V.  Using  a  discrete  version  of 
equation  (55),  our  goal  is  to  establish  the  integral  representations  in  (28)  for  the  functionals 
S*fc  =  e(5jk  +  (XjiXk)  1,2)  and  A*fc  =  (iMxj,Xk)  1,2  in  (48),  involving  a  discrete  spectral 
measure  which  is  analogous  to  the  measure  in  equation  (34). 

By  our  discussion  in  Section  IV  A,  the  matrix  representation  of  (55)  is  given  by 

Bzn  —  AnCzn,  B  =  —/VIHV,  C  =  VTV.  (56) 

The  first  formula  in  equation  (56)  is  a  generalized  eigenvalue  problem  [40]  associated  with 
the  pencil  B  —  AC,  where  B  and  C  are  Hermitian  and  real-symmetric  matrices,  respectively, 
of  size  K,  say.  The  Xn  and  zn,  n  =  1, . . . ,  K,  are  known  as  generalized  eigenvalues  and 
eigenvectors,  respectively.  The  matrix  C  =  V7  V  is  clearly  positive  semi-definite.  In  this 
section,  we  will  assume  that  C  is  positive  definite,  hence  invertible.  We  will  discuss  the  case 
where  C  is  positive  semi-definite  in  Section  VIII. 

Since  the  matrices  B  and  C  are  symmetric  and  C  is  positive  definite,  the  generalized 
eigenvalue  problem  has  properties  which  are  similar  to  the  properties  of  the  standard  sym¬ 
metric  eigenvalue  problem  [40].  In  particular,  the  generalized  eigenvalues  An  are  all  real,  the 
generalized  eigenvectors  zn  form  a  basis  for  CA ,  and  the  zn  are  orthonormal  in  the  following 

sense  z,[Czm  =  8nrn ,  n,m  —  1, _ ,  K,  [40].  Since  C  =  V7  V  is  real- valued,  this  is  equivalent 

to  the  Sobolev-type  orthogonality  condition 

V zm  =  8nm.  (57) 

In  other  words,  the  generalized  eigenvectors  zn  are  orthonormal  with  respect  to  the  “discrete 
inner-product”  (-,  -)i,2  defined  by  (£,C)i,2  =  V£-V£,  for  £  G  CA.  Denoting  by  Z  the 
matrix  with  columns  consisting  of  the  generalized  eigenvectors  zn ,  equation  (57)  is  seen  to 
be  equivalent  to  [VZ]1[VZ]  =  I,  or  Z^CZ  =  I.  A  key  feature  of  the  generalized  eigenvalue 
problem  is  that  the  matrix  Z  simultaneously  diagonalizes  B  and  C.  Specifically,  if  A  is  the 
diagonal  matrix  whose  elements  on  the  main  diagonal  are  the  generalized  eigenvalues  An, 
then  [40] 


ZfBZ  =  A, 


ZfCZ  =  I. 


(58) 
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Since  the  zn,  n  —  1, . . . ,  K,  form  a  basis  for  CA  and  satisfy  (57),  for  all  £  G  CA,  we 
have  that  $,  =  ^n(Vzn*V£)zn  =  Y2n(zn[^zn]  ^V)£,  which  implies  the  following  analogue 
of  equation  (33) 

K 

Y,  Qn  =  I,  Qn  =  Zn[Vz„]  fV,  QlQm  =  Ql  (59) 

71=1 

where  the  matrices  Qn,  n  =  1, . . . ,  K,  are  self-adjoint  with  respect  to  the  discrete  inner- 
product  (•,  -)i,2  defined  above,  i.e.,  (Qn£,  C)i,2  =  (£,  QnC)i,2  for  all  £,  £  G  CA. 

We  now  use  equation  (59)  to  prove  the  spectral  theorem  in  (22)  for  this  generalized  eigen¬ 
value  problem  setting.  From  B zn  =  AnC zn  and  equation  (59)  we  have  that  B  Qn  =  AnC  Qn 
which  is  equivalent  to  C_1BQn  =  A nQn,  since  the  matrix  C  is  invertible.  As  in  the  discussion 
following  equation  (33),  the  mutual  orthogonality  of  the  projection  matrices  Qn  implies  that 
/(C_1B)  =  Y^n /(An)Qn  for  any  polynomial  /  :  1  h  C.  From  the  mutual  orthogonality  of 
the  projection  matrices  Qn  and  their  symmetry  with  respect  to  the  discrete  inner-product 
(•,  it  follows  that,  for  all  £,  £  G  CA  and  complex-valued  polynomials  /(A)  and  h( A),  the 
bilinear  functional  (/(C_1B)^,  h(C_1B)^)ii2  has  the  integral  representation  in  (22),  with  M 
substituted  by  C_1B.  Moreover,  the  complex-valued  function  fJ,^(X)  =  (Q(A)£,£)  in  (22) 
is  now  given  by  /i^(A)  =  (Q(A)£,  C)i,2,  where  the  associated  matrix  representation  Q(A) 
of  the  projection  operator  Q(A)  and  the  discrete  spectral  measure  d/i^(A)  are  given  by  the 
following  analogue  of  equation  (34) 

S(A)  =  £  0(A-A„)Q„,  dftc(A)=  J2  <«A»(dA)[VQ„«-VC]>.  (60) 

n:  An<A  n:  An<A 

We  now  demonstrate  that  discrete  representations  of  the  functional  formulas  for  S*k 
and  A*k  in  equation  (48)  yield  their  associated  integral  representations  in  (28),  involving 
the  discrete  spectral  measure  in  (60).  From  A  =  A_1(V-HV)  and  (56),  the  matrix 
representation  of  the  functional  formulas  S*fc  =  s(Sjk  +  (xj,  Xk)  1,2)  and  A*k  =  (Axj,  Xk)  1,2  in 
equation  (48)  are  given  by 

s*„  =  Efe  +  (VXj-VXk))  A*t  =  (61) 

Moreover,  the  matrix  representation  of  the  cell  problem  eA \j  +  V-HVyq  =  —Uj,  in  (4),  is 
given  by 


(cC  +  iB)Xj  =  Uj, 


(62) 
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where  Uj  is  the  discrete,  vector  representation  of  the  jth  component  of  the  fluid  velocity 
held  Uj,  for  example.  The  matrix  Z  is  invertible,  as  the  generalized  eigenvectors  zn  form  a 
basis  for  CN.  Consequently,  by  equation  (58)  we  have  that  B  =  Z_lAZ_1,  C  =  Z  ^Z -1.  It 
now  follows  from  equation  (62)  that  Z~\e\  +  iA)Z~1Xj  =  Uj,  or  equivalently 

Xj  —  Z(s\  +  (63) 

Substituting  the  resolvent  formula  for  Xj  in  (63)  into  equation  (61)  yields  the  following 
formula  that  is  a  direct  analogue  of  equation  (36) 

S*jk  =  e{5jk  +  ((el  +  ?A)_1ZtUj*(el  +  *A)_1ZTwfc))  (64) 

A*k  =  (*A(el  +  ?A)_1Z1'uj*(el  +  /A)  'Z 

where  we  have  used  that  [VZ]t  =  [VZ]_1.  The  quadratic  form  Z^uj-Z^uk  arising  in  (64) 
can  be  written  in  terms  of  the  projection  matrices  Qn  defined  in  (59)  as  follows.  Analogous 
to  equation  (37),  we  have 

N  _  N 

Z]Uj‘Z]Uk  =  22  ( ziuj)(z]iUk )  =  22  ZnZ]nUj-Uk.  (65) 

n=  1  n=  1 

We  now  demonstrate  that  znzJ2Uj-uk  =  V QnflyV<7fe,  where  gJ  =  (VTV)_1Wj, 

znz\uj-uk  =  znz\ [VT V]^ •  [VTV]grfc  =  [Vzn][Vzn]  ]Vgj-Vgk  =  VQ„^*V^fc.  (66) 

It  follows  from  equations  (65)  and  (66)  that  the  functional  formulas  for  S*k  and  A*k  in  (64) 
have  the  Stieltjes  integral  representations  in  equation  (26),  involving  the  discrete  spectral 
measure  / ijk  in  (60)  with  £  =  <7  ■  =  (V7  V)-1Mj  and  C,  =  9k  =  (VTV)-1ttfc- 

Exactly  as  in  Section  IV  A,  we  may  rewrite  these  integral  representations  for  S*fc  and 
A*k  in  equation  (26),  involving  the  complex  measure  gjk,  as  the  integral  representations  in 
equation  (28)  involving  the  real  signed  measures  Re  gjk  and  Im  /i^.  The  associated  real¬ 
valued  functions  Re  pjk( A)  and  Im  /j,jk( A)  are  given  by  equation  (38)  with  [(Qn  +  Qn)gj-gk] 
substituted  by  [V(Qn  +  Qn)9jmVgk].  Furthermore,  due  to  the  generalized  eigenvalues  and 
eigenvectors  of  the  antisymmetric  matrix  ?.M  =  V7  HV  coming  in  complex  conjugate  pairs, 
the  discrete  measure  /.ijk  depends  only  on  a  subset  of  the  index  set  of  the  sums  in  (65)  and 
satisfy  equation  (39),  with  [Qngj-gk  ]  substituted  by  [VQn<jyVflffc  ]. 
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VII.  DISCRETE  EQUIVALENCE  OF  THE  EFFECTIVE  PARAMETER 

PROBLEMS 

In  Section  III,  we  formulated  an  effective  parameter  problem  for  the  effective  diffusivity 
tensor  D*  associated  with  an  incompressible  fluid  velocity  held.  The  discrete  version  of 
this  effective  parameter  problem  was  formulated  in  Section  IV.  An  alternate  approach  to 
the  effective  parameter  problem  was  formulated  in  Section  V,  and  its  discrete  version  was 
formulated  in  Section  VI.  In  this  section,  we  demonstrate  that  the  discrete  versions  of  these 
effective  parameter  problems  yield  equivalent  spectral  representations  of  D*  when  the  matrix 
V  is  of  full  rank,  so  that  the  matrix  Laplacian  is  invertible. 

Let  V  =  UEV7  be  the  singular  value  decomposition  (SVD)  of  the  matrix  V  of  size  N  x  K, 
say,  where  N  >  K.  Here,  E  =  diag(<7i, . . . ,  <Jk),  where  0  <  o\  <  ■  ■  ■  <  ax,  and  the  matrices 
U  and  V  are  of  size  N  x  K  and  K  x  K,  respectively,  and  satisfy  [13] 

UTU  =  I,  VTV  =  VVT  =  I,  (67) 

where  I  is  the  K  x  K  identity  matrix.  The  columns  of  U  are  called  left  singular  vectors,  the 
columns  of  V  are  called  right  singular  vectors,  and  the  a*  are  called  singular  values. 

It  follows  from  V  =  UEV7  and  equation  (67)  that  the  spectral  decomposition  of  the 
negative  matrix  Laplacian  VTV  is  given  by  VTV  =  VE2V7  [13].  We  assume  that  V  is  of 
full  rank  so  that  cq  >  0  for  all  %  —  1, . . . ,  K.  This  implies  that  LV1  exists  so  that  the  matrix 
Laplacian  is  invertible.  In  this  case,  it  follows  from  V  =  UEV7  and  equation  (67)  that  the 
projection  matrix  T  =  V(VrV)_1Vi  is  given  by 

T  =  UUT,  (68) 

which  is  a  N  x  N  symmetric  projection  matrix  satisfying  T2  =  T  and  TV  =  V  (see  equa¬ 
tion  (67)).  A  key  property  of  the  SVD  of  the  full  rank  matrix  V  is  that  its  range  is  spanned 
by  the  columns  of  U  [13],  hence  T  =  UU7  projects  subspaces  of  onto  the  range  of  V. 

From  equations  (67)  and  (68),  we  can  write  the  eigenvalue  problem  — ?rHTiun  =  A nwn 
discussed  in  Section  IV  A  as 

[— *UtHU][Utiu„]  =  An[UTiun].  (69) 

Now  consider  the  generalized  eigenvalue  problem  —  ?VTHVzn  =  «nVJ  Vz„  discussed  in 
Section  VI  and  recall  that  V  =  UEV7  and  V7V  =  VS2V7 .  Since  E  is  invertible,  by 
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equation  (67),  we  can  write  this  generalized  eigenvalue  problem  as  the  following  standard 
eigenvalue  problem 

HUrHU][E  VTzn]  =  an[YMTzn\.  (70) 

Comparing  the  formulas  in  equations  (69)  and  (70)  indicates  that  spectrum  associated  with 
each  of  these  eigenvalue  problems  is  identical,  an  =  Xn,  and  that  the  eigenvectors  are 
related  by  liTwn  =  EVT zn.  Since  T  is  a  projection  matrix,  T2  =  T,  the  eigenvalue  problem 
FHTwn  =  i\nwn  can  be  written  as  rHr[ru;n]  =  *An[riyn],  which  implies  that  Vwn  =  wn. 
Consequently,  applying  the  matrix  U  to  both  sides  of  the  formula  U7  wn  =  EV7  zn  and 
recalling  that  T  =  UU7  and  V  =  UEVr  we  have 

wn  =  Vzn.  (71) 

In  the  following  lemma  we  make  precise  the  correspondence  between  the  standard  eigen¬ 
value  problem  —  ?rHrtun  =  A nwn  and  the  generalized  eigenvalue  problem  — ?VTHVzn  = 
a„VTVz„,  as  well  as  the  associated  spectral  measures  in  equations  (34)  and  (60),  respec¬ 
tively. 

Lemma  1  Consider  the  standard  eigenvalue  problem  and  the  generalized  eigenvalue  problem 
given,  respectively,  in  equations  (72)  and  (73)  below 

—zTHTwn  =  X  nwn,  (72) 

— ?VTHVzn  =  XnVTVzn.  (73) 

Let  V  =  UEV7  be  the  SVD  of  the  matrix  V,  which  we  assume  to  be  of  full  rank.  Then 
equation  (72)  implies  and  is  implied  by  equation  (73),  with  wn  and  zn  related  as  in  equa¬ 
tion  (71) .  This  implies  that  the  spectrum  associated  with  each  of  these  eigenvalue  problems 
is  identical.  Moreover,  the  spectral  weights  in  equations  (37)  and  (66)  are  identical,  i.e., 

QnrHerrHefc  =  VQn[VrV]- VV[VTV]“W  (74) 

This,  in  turn,  implies  that  the  associated  spectral  measures  in  equations  (34)  and  (60)  are 
identical  for  all  £  G  C,v. 

Proof  of  Lemma  1 

Recall  that  V  =  UEVr,  VTV  =  VE2Vr,  and  V  =  UU7  ,  where  E  is  invertible,  and  the  matri¬ 
ces  V  and  U  satisfy  equation  (67).  First  consider  equation  (72)  written  as  in  equation  (69), 
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[ — zllT  HU][UTtun]  =  An[UTtun].  Since  the  matrix  E  is  invertible  and  V7  V  =  I,  we  can  rewrite 
equation  (69)  as 

VE[-2UtHU](EVt)(VE-1)[Ut^„]  =  An(VE2VT)  (VE_1)  [\JTwn\ ,  (75) 

which  is  precisely  equation  (73)  written  in  terms  of  V  =  UEV7  with  zn  =  VE_1UTmn.  This 
formula  for  zn,  equation  (67),  and  the  formula  Twn  =  wn  above  equation  (71)  imply  that 
wn  =  UE MT zn  =  Vzn.  Now  consider  equation  (73)  written  as  in  (70),  [— ?U7  HU]  [EV7  zn\  = 
\n[YMTzn}.  Since  UTU  =  I,  we  can  rewrite  equation  (70)  as 

U[— ?UTHU](UTU)[EVTz:n]  =  AnU[E  VTzn],  (76) 

which  is  precisely  (72)  written  in  terms  of  T  =  UU7  with  wn  =  UE \/T zn  =  Vz„. 

We  now  establish  equation  (74).  From  the  formula  u  =  V-H  in  (8),  we  have  that 
Uj  =  V-H ej.  Since  V  =  UEVT,  the  discrete  version  of  this  formula  is  given  by 

uj  =  -VTH  ej  =  -VEUTHei.  (77) 

From  V  =  UEVT  and  (VTV)-7  =  VE“2VT  we  have  V(VTV)"1  =  UE-W7  Consequently, 
T  =  UUr,  and  equations  (67)  and  (77),  yield  — FH ej  =  V(VTV)_1Wj.  Equation  (74)  now 
follows  from  the  formula  wn  =  Vzn  in  (71) 

wnwfn v  (' VTV) ■ V  (' VTV) ^Uk  =  [(' VT' V)  ■ -1' VTwn]  [( VTV)  -1  VTru„]  -uk 

=  [(VTV)-1VTVzn][(VTV)-1VTVzn]t'ui-'ufc 
=  Znziuj-Uk,  (78) 

where  we  have  used  that  the  inverse  of  a  symmetric  matrix  is  also  symmetric  [22],  The 
equivalence  of  equations  (74)  and  (78)  now  follows  from  equations  (33),  (59),  and  (66).  This 
concludes  our  proof  of  Lemma  1  □ . 

We  conclude  this  section  with  a  discussion  regarding  numerical  computations  of  the  ef¬ 
fective  diffusivity  D*.  The  approach  discussed  in  this  section  and  the  projection  method 
discussed  in  Section  IV  B  combine  the  computational  advantages  of  the  methods  discussed 
in  Sections  IV  A  and  VI.  More  specifically,  in  the  full  rank  setting,  the  spectral  measure 
underlying  the  discrete  integral  representation  for  D*  was  calculated  in  Section  IV  A  in 
terms  of  the  standard  eigenvalue  problem  — *rHrmm  =  A nwn,  where  the  matrix  — ^THT 
is  of  size  N  x  N.  In  Section  VI,  D*  was  calculated  in  terms  of  the  generalized  eigenvalue 
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problem  — ?V7  H V zn  =  A nV7  Vzn,  involving  the  K  x  A'  matrices  —  iSJ1  H  V  and  V7  V.  Since 
VT  =  (Vf, . . . ,  V  J)  is  of  size  K  x  N  we  have  that  A  =  N/d.  However,  the  generalized 
eigenvalue  problem  is  more  computationally  intensive  than  the  standard  eigenvalue  prob¬ 
lem  [40].  For  the  case  of  randomly  perturbed  flows,  these  are  not  efficient  ways  to  compute 
spectral  statistics  for  D*. 

The  projection  method  developed  in  Section  IV  B  demonstrates  that,  by  first  computing 
the  standard  eigenvalue  decomposition  of  the  non-random  matrix  T,  the  spectral  statistics 
of  the  eigenvalue  problem  —  ?rHTiun  =  Xnwn  can  then  be  obtained  by  repeatedly  comput¬ 
ing  the  standard  eigenvalue  decomposition  of  smaller  matrices.  They  are  of  size  A  x  A 
by  equations  (40),  (41),  and  (68).  We  stress  that  A  and  N\  both  denote  the  rank  of  the 
matrix  T  in  this  section  and  in  Section  IV  B,  respectively,  i.e.,  A  =  N\ .  Note  that  com¬ 
puting  the  matrix  T  =  V(VTV)_1V7  involves  numerically  solving  N  linear  systems  of  size 
A  x  A.  Alternatively,  the  proof  of  Lemma  1  illustrates  that  by  first  computing  the  SVD  of 
the  matrix  gradient,  V  =  UEVT,  the  spectral  statistics  of  the  generalized  eigenvalue  prob¬ 
lem  — ?,V7  H W zn  =  AnV7  V zn  can  then  be  obtained  by  repeatedly  computing  the  standard 
eigenvalue  decomposition  of  the  matrix  —  ?U7  HU  which  is  of  size  A  x  K .  When  N  is  large, 
these  equivalent  methods  are  more  numerically  efficient  than  the  other  approaches  discussed 
in  Sections  IV  A  and  VI. 

In  Section  VIII,  we  generalize  Lemma  1  to  the  case  where  V  has  rank  K\  with  K\  <  K , 
demonstrating  that  the  two  formulations  yield  equivalent  spectral  representations  of  the 
effective  diffusivity  tensor  D*  in  the  this  rank  deficient  setting.  Moreover,  we  demonstrate 
that,  by  first  computing  the  SVD  of  the  matrix  gradient,  D*  can  be  computed  via  a  standard 
eigenvalue  problem  for  matrices  of  size  K\  x  K\.  Consequently,  the  rank  deficiency  of  the 
problem  actually  increases  the  numerical  efficiency  of  computations. 


VIII.  RANK  DEFICIENCY  AND  A  UNIFYING  STANDARD  EIGENVALUE 

PROBLEM 

In  Sections  IV  and  VI  we  provided  two  discrete,  matrix  formulations  of  the  effective 
parameter  problem  for  the  effective  diffusivity  tensor  D*.  These  two  formulations  assume 
that  the  N  x  K  matrix  V  is  of  full  rank  K  so  that  the  negative  matrix  Laplacian  V7  V  is 
invertible.  Lemma  1  shows  that,  given  this  condition,  the  two  formulations  yield  equivalent 
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spectral  representations  of  D*.  In  this  section,  we  generalize  Lemma  1  to  the  case  where  V 
has  rank  K\  with  K\  <  K,  again  demonstrating  that  the  two  formulations  are  equivalent  in 
this  rank  deficient  setting.  This  framework  is  used  in  Section  IX  to  compute  D*  for  periodic 
flows,  for  which  the  matrix  V  with  periodic  boundary  conditions  is  rank  deficient. 

Consider  the  cell  problem  in  (4)  written,  via  (8)  and  [V*H]*V</?  =  V*[HV</?],  as 

V-HVxj  +  sA  Xj  =  -uj.  (79) 

Discretizing  equation  (79)  yields  (see  Section  IV  A  for  details) 

VTH  VXj  +  =  uj ,  (80) 

where  Uj  is  the  discrete,  vector  representation  of  the  jth  component  of  the  fluid  velocity  field 
Uj,  and  similarly  for  Xj-  Substituting  the  formula  for  Uj  in  (80)  into  the  discrete  version 
D*k  =  eSjk  +  (uj'Xk)  °f  equation  (3)  yields 

D-jh  =  S-k  +  A-t,  S'jk  =  E(6jt  +  (VTVXrx*)),  A*t  =  (VTH  VXj-Xt),  (81) 

where,  as  before,  S*kj  =  S*k  and  A *k-  =  —  A*fc. 

We  first  work  with  equation  (80)  directly  and  develop  a  mathematical  framework  which 
parallels  the  framework  of  Section  VI.  We  then  transform  equation  (80)  into  a  discrete  ana¬ 
logue  of  equation  (18)  written  as  (el  +  THT)  VXj  —  —  THej,  with  a  suitable  generalization  of 
the  formula  for  the  matrix  T  in  (68) ,  and  develop  a  mathematical  framework  which  parallels 
the  framework  of  Section  IV.  We  then  generalize  Lemma  1  of  Section  VII,  establishing  the 
equivalence  of  these  two  formulations  for  the  rank  deficient  setting. 

Let  V  =  USVT  be  the  SVD  of  the  N  x  K  matrix  V  discussed  in  Section  VII.  We  assume 
that  V  is  rank  deficient  so  that  VTV  =  VX2V7  is  singular,  with  K\  non-zero  eigenvalues 
and  Kq  =  K  —  K\  zero  eigenvalues,  and  write 

U  =  [U0  Ux],  X  =  °00  °01  ,  V  =  [V0  Vx].  (82) 

Oio  Si 

Here,  0 ab  are  matrices  of  zeros  of  size  Ka  x  Kb,  a,b  =  0, 1,  Ua  is  N  x  Ka ,  Va  is  K  x  Ka , 
and  Si  is  a  K\  x  K\  diagonal,  invertible  matrix.  By  equation  (67)  the  matrices  Ui  and  Vi 
satisfy 


ufu,  =  I,.  Vjv,  =  I,, 


(83) 
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where  h  is  the  K\  x  K\  identity  matrix.  An  important  property  of  the  SVD  of  the  matrix 
V  is  that  its  null  space  is  spanned  by  the  columns  of  Vo  and  its  range  is  spanned  by  the 
columns  of  Ui  [13]. 

Due  to  the  blocks  of  zeros  of  E  in  (82),  the  matrix  elements  of  V  and  V7  V  do  not  depend 
on  Uo  nor  V0  and  can  be  written  as  V  =  UxExV7  and  V7  V  =  VxE7V7.  Consequently, 
equation  (80)  can  be  rewritten  as 

[V1E1][Uj’HU1][E1Vf]x3  +  eVjEXXi  =  nr  (84) 

Since  the  K\  x  I\\  matrix  U  ( HU  i  is  antisymmetric,  it  has  the  spectral  decomposition 
UfHUi  =  zRiAiR},  where  Ai  is  a  diagonal,  real-valued  matrix  and  Ri  is  a  unitary  matrix, 
R{Ri  =  RiR|  =  1 1 .  It  follows  that 

VrHV  =  *[V1E1R1]A1[V1E1R1]t,  VTV  =  [VxEiRiHVxExRx]*.  (85) 

Consequently,  equation  (84)  can  be  rewritten  as  [ViEiRi](di  -MAi) [ViEiRi^Xj  =  uj-  This 
formula  and  equation  (83)  together  imply  the  following  analogue  of  equation  (63) 

VfXi  =  VfZxOHx  +iAl)~1Z{uj,  Zx  =  VxE^Rx.  (86) 

Substituting  equation  (86)  into  the  functionals  {VT\/Xj'Xk)  and  (VT  HVxpXfc)  hr  equa¬ 
tion  (81)  yields  the  following  analogue  of  (64)  (see  Appendix  B  for  details) 

S*fc  =  £(Sjk  +  ((eli  +  *Ax)^1Zjwi-(£li  +  zAx)_1Z  |  itfc)),  (87) 

A*fc  =  (*Ax(elx  +  *Ax)_1Zjwi-(eli  +  iA1ylZ\uk). 

The  quadratic  form  Z\uj-Z\uk  arising  in  (87)  yields  the  following  analogue  of  (65) 

z\urz\uk  =  ([z^u^dz^Uk)  =  [Zn]  ^j-Uk,  z\  =  VjE 7Y7 ,  (88) 

n= 1  n=  1 

where  r7,  n  =  1, _ ,Ki,  are  the  orthonormal  eigenvectors  of  the  matrix  UfHUx  which 

comprise  the  columns  of  Rx.  From  V  =  UxExVf  and  equations  (83)  and  (88)  we  have 
that  Vz*  =  Ux r7.  The  orthogonality  condition  r^-r^n  =  8nm  and  equation  (83)  then 
imply  that  the  vectors  z\  satisfy  the  Sobolev-type  orthogonality  condition  in  equation  (57), 
Vz*-VZjn  =  Snm .  Moreover,  since  =  Jx,  as  the  vectors  r]n  form  an  orthonormal 

basis  for  CKl,  we  also  have  the  following  generalization  of  equation  (59) 

Ki 

^Ci  =  v,vf,  Qi  =  zi\Vzi  ]fV,  Q}Q'm  =  Q}Slm,  (89) 

n=  1 
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where  the  matrices  Q*,  n  —  1, . . .  ,  are  self-adjoint  with  respect  to  the  discrete  inner- 

product  (•,  * )  1,2  defined  by  (£,  C)i,2  =  (V£-VC),  i.e.,  (Q„£,  C)i,2  =  (£,  QiC)i,2  for  £,  C  e  CKl. 
It  follows  from  equations  (87)  and  (88)  that 


Ki 


S-Je  -  Sjt  =  £ 


Re 


u 


:)([*!]'«*)  ] 


n=  1 


+  [^n] 


A*  _  'W' 


u 


0Mf«*)] 


n=l 


+  WJ" 


(90) 


where  A*,  n  —  1, . . . ,  K\ ,  are  the  eigenvalues  of  the  matrix  — ?UfHUi  corresponding  to  the 
eigenvectors  r\.  Here,  as  in  equation  (28),  we  have  used  the  fact  that  the  matrices  V  and 
H,  as  well  as  the  vectors  Xi  and  uji  and  the  molecular  diffusivity  £  are  real  valued,  so  that 
(VXj-VXfc)  =  (VXfc'VXj)  and  (VTHVxi-%fc)  =  (xfc*VrHVXj).  One  to  the  eigenvalues 
and  eigenvectors  of  the  antisymmetric  matrix  Uf  HUi  coming  in  complex  conjugate  pairs,  as 
in  equations  (38)  and  (39),  the  sums  in  (90)  depend  only  on  a  subset  of  the  index  set. 

We  now  argue  that  the  mathematical  framework  developed  in  equations  (80)-(90)  gen¬ 
eralizes  the  full  rank  case  in  Section  VI.  Indeed,  in  the  full  rank  setting,  the  matrix  Vi  =  V 
is  orthogonal,  Hi  =  E  is  invertible,  and  Ri  =  R  is  orthogonal,  so  that  the  matrix  Zi  =  Z 
defined  in  (86)  is  given  by 


Z  =  VE_1R, 


(91) 


and  is  invertible  with  Z_1  =  R^EV7  and  satisfies  Z  XZ  =  ZZ  1  =  I.  Consequently,  equa¬ 
tion  (85)  implies  that  equation  (58)  is  satisfied  with  Ai  =  A.  In  this  case,  equations  (63)-(65) 
are  identical  to  equations  (86)-(88) ,  respectively. 

We  now  generalize  the  mathematical  framework  developed  in  Section  IV  to  the  case  that 
the  matrix  V  is  rank  deficient.  Using  equation  (83)  and  the  invertibility  of  the  matrix  Ei, 
we  can  rewrite  equation  (84)  as 

U^UfHU^UfU^EW^x,  +eU1X1V[xj  =  UiE^Vfu,,  (92) 

Substituting  the  formula  Uj  =  — VrHej  of  (77)  into  equation  (92)  and  using  V  =  UiEiVf 
yields 


(e\  +  r1Hr1)vxj  =  g],  r\  =  ihuf,  ^-ihHe,,  (93) 

which  analogous  to  equation  (18).  As  in  Section  VII,  the  matrix  Y \  =  IRUf  projects 
subspaces  of  WN  onto  the  range  of  V.  Since  the  matrix  IAHTi  is  antisymmetric,  it  has 
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the  spectral  decomposition  Id  Hid  =  AA/iAWj,  where  Wi  is  a  unitary  matrix  WjW,  = 
WiW{  =  I.  This  and  equation  (93)  yield  the  resolvent  formula  for  VXj  hr  (35),  with 
corresponding  notational  changes.  In  equation  (81),  write  (VT\/Xj‘Xk)  =  (^Xj'^Xk)  and 
(VTHVXj-Xfc)  =  (TxH  TiVx/VXfc),  where  the  second  formula  follows  from  UfUi  =  h 
and  V  =  U1E1V1  which  imply  that  TiV  =  V.  Exactly  as  in  Section  IV  A,  this  leads  to 
equations  (36)  and  (37),  with  corresponding  notational  changes.  This,  in  turn,  leads  to  the 
integral  representations  for  S*k  and  A*fc  in  equation  (28) ,  involving  a  discrete  spectral  measure 
Hjk  associated  with  the  function  Pjfc(A)  =  (QX(A )g)-g)),  defined  by  (34)  with  Q(A)  and  Qn 
substituted  by  QX(A)  and  Q7,  respectively,  where  Q7  =  and  wln ,  n  —  1, . . .  ,Ad, 

are  the  eigenvectors  of  the  matrix  TiHTi  which  comprise  the  columns  of  Wi.  In  the  case 
that  the  matrix  V  is  of  full  rank,  we  have  Ui  =  U  hence  Wi  =  W.  This  establishes  that  the 
mathematical  framework  discussed  in  this  paragraph  reduces  to  the  mathematical  framework 
in  Section  IV  for  the  full  rank  setting  and  therefore  is  a  generalization. 

We  now  employ  the  projection  method  developed  in  Section  IV  B  to  generalize  the  math¬ 
ematical  framework  in  Section  VII  to  the  rank  deficient  setting,  establishing  the  equivalence 
of  the  two  approaches  that  follow  from  equations  (84)  and  (93).  We  stress  that  the  matrix 
Ti  defined  in  this  section  is  identical  to  the  matrix  T  defined  in  Section  IV  B,  Ti  =  T,  which 
were  given  different  notations  to  clarify  our  discussion  here. 

Since  Ti  =  U  i  Uf  is  a  projection  matrix,  the  projection  method  discussed  in  Section  IV  B 
can  be  directly  applied  to  equation  (93).  However,  the  mathematical  framework  developed 
here  provides  deeper  insight  into  equation  (41)  of  the  method.  In  particular,  in  Section  IV  B 
we  wrote  T  =  PGP7,  where  P  is  an  orthogonal  matrix  consisting  of  the  eigenvectors  of  T 
and  the  matrix  G  is  defined  in  equation  (41).  Moreover,  we  wrote  P  =  [Po  Pi],  where  the 
columns  of  the  matrices  Po  and  Pi  are  orthonormal  eigenvectors  that  span  the  null  space 
and  range  of  T,  respectively.  Since  the  eigenvalues  yn  associated  with  the  eigenvectors  in  the 
matrix  Pi  satisfy  =  1,  any  linear  combination  of  the  corresponding  eigenvectors  is  also 
an  eigenvector  of  T  with  eigenvalue  yn  =  1.  Therefore,  since  the  orthonormal  columns  of  the 
matrix  Ui  span  the  range  of  Ti  =  UiUf,  without  loss  of  generality,  we  may  take  Pi  =  Ui  so 
that  P=[P0Ui],  Consequently,  we  can  rewrite  equation  (42)  as 

UfHUi^.RiiAnR^.  (94) 

From  equation  (94)  and  the  comment  after  equation  (84),  we  have  that  Rn  =  Ri  and 
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An  =  Ai.  This  and  equation  (43)  establishes  that  the  spectra  associated  with  each  of  the 
two  approaches  are  identical.  We  now  establish  that  the  spectral  weights  associated  with 
both  approaches  are  also  identical.  From  the  formula  P  =  [Po  Ui]  and  equation  (43)  with 
Rn  =  Ri  and  W  redefined  as  Wi,  we  have  that  Wi  =  PR  =  [Po  Ui Ri] .  Consequently, 
from  V  =  UiEiVf,  equation  (83),  and  the  formula  Zi  =  ViE^Ri  in  (86),  we  have  that 
VZi  =  UiRi,  which  implies  the  following  generalization  of  equation  (71) 

Wi  =  [P0  VZi],  VZi  =  UiRi.  (95) 

It  now  follows  from  TiPo  =  0,  TiV  =  V,  and  the  formula  Uj  =  —  V7  He.,  in  (77)  that 

WlriHej-WjriHefc  =  [VZ^He^VZ^Hefc  =  Z[urZ\uk.  (96) 

This  establishes  that  the  spectral  weights  associated  with  both  approaches  are  identical  and, 
in  turn,  establishes  the  equivalence  of  the  two  approaches  following  from  equations  (84) 
and  (93).  In  Section  IX  we  will  use  equation  (90)  to  compute  the  effective  diffusivity  ten¬ 
sor  D*  for  various  periodic  flows  and  relate  spectral  characteristics  to  flow  geometry  and 
transport  properties. 


IX.  COMPUTATIONS  OF  THE  EFFECTIVE  DIFFUSIVITY  TENSOR 

In  this  section,  we  employ  the  mathematical  framework  developed  in  Section  VIII  to  pro¬ 
vide  rigorous  computations  of  the  effective  diffusivity  tensor  D*  for  various  model  periodic 
flows.  In  particular,  we  employ  equation  (90)  to  compute  the  symmetric  S*  and  antisymmet¬ 
ric  A*  parts  of  D*.  As  a  benchmark  test  we  compute  D*  for  shear  flow,  for  which  the  spectral 
measure  is  known  [3] .  We  also  consider  a  fluid  velocity  field  that  has  a  free  parameter.  As 
the  parameter  varies,  the  flow  transitions  from  cellular  flow  to  shear  flow  in  the  diagonal  x-y 
direction.  This  gives  rise  to  fascinating  transitional  behavior  in  the  spectral  measure,  which 
governs  transitional  behavior  in  the  effective  diffusivity  tensor  D*.  For  the  sake  of  brevity, 
we  will  focus  our  attention  on  the  £  behavior  of  the  components  S*fc,  j,  k  —  lt , . . ,  d,  of  S*. 
Moreover,  for  computational  simplicity,  we  have  restricted  our  computations  to  dimension 

d  =  2. 
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A.  Numerical  Methods 


By  equation  (8),  the  time-independent  fluid  velocity  field  u  =  u(x)  is  given  in  terms 
of  an  antisymmetric  matrix  H  =  H(.x),  u  =  V-H.  For  dimension  (1  =  2.  the  matrix  H  is 
determined  by  a  stream  function  'I'  =  T(:e), 


0  T 
-T  0 


(97) 


yielding  u  =  [— <9X'F],  where  dx  denotes  partial  differentiation  in  the  variable  x,  for 
example.  In  this  section  we  consider  two  flows  with  free  parameters  which  transition  from 
cellular  flow  to  shear  flow  as  parameters  vary.  In  particular,  we  consider  BC-flow  [9]  and 
“cat’s  eye”  flow  [15],  which  are  defined  by  the  following  stream  functions  'F bc  and  'F ce , 
respectively, 


^bc(x)  —  Bsinx  —  Csiny,  ^ce (%)  =  sin x  sin y  +  A  cos x  cosy,  (98) 

where  we  have  denoted  x  =  ( x,y ).  The  corresponding  fluid  velocity  fields  are 

ubc{x)  =  (C  cos  y,  B  cos  a;),  (99) 

uce  (*)  =  (—  sin  x  cos  y  +  A  cos  x  sin  y,  cos  x  sin  y  —  A  sin  x  cos  y) . 


The  flow  geometry  of  these  fluid  velocity  fields  transition  from  shear  to  cellular  flow  structure 
as  the  system  parameters  A,  B,C  €  [0, 1]  vary. 

Streamlines  of  a  2D  flow  are  level  sets  of  the  stream  function  T,  which  define  a  family  of 
curves  that  are  instantaneously  tangent  to  the  fluid  velocity  field  u,  since  u  =  [— dy^,  <9XVF] 
implies  that  w-VT  =  0.  In  Fig.  1,  we  display  streamlines  of  the  stream  functions  in 
equation  (98)  for  various  parameter  values.  The  streamlines  for  cat’s  eye  flow  are  sym¬ 
metric  about  the  line  y  =  x,  which  follows  from  the  symmetry  of  the  stream  function 
VI ’cE(x,y)  =  \F cE{y,x )•  The  stream  function  for  BC-flow  has  a  more  complicated  sym¬ 
metry  ty(x,y,  B,C)  =  —  \F(?/,  x;  C,  B).  This  symmetry  indicates  that  if  the  values  of  B 
and  C  are  interchanged,  B  < — »  C ,  then  the  original  flow  is  recovered  from  a  90°  rotation 
(x  — >  y,  y  — »  —x)  followed  by  a  reflection  about  the  x-axis  (y  — >  —y).  Consequently,  flows 
elongated  in  the  ^/-direction  become  flows  elongated  in  the  x-direction  under  the  interchange 

B< — >  C. 
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Cat's  Eye  Flow 


FIG.  1.  Transitions  in  flow  structure.  The  streamlines  for  BC-flow  and  cat’s  eye  flow  are  displayed 
for  various  parameter  values,  transitioning  from  shear  to  cell  flow  structure.  From  left  to  right 
and  top  to  bottom,  the  parameter  values  associated  with  BC-flow  are  B  =  1  fixed  and  C  = 
0,  0.01,  0.1,  0.3,  0.5  and  1,  while  those  for  cat’s  eye  flow  are  ^4  =  0,  0.1,  0.3,  0.5,  0.7,  and  1. 

In  equation  (6)  and  the  paragraph  therein  we  discussed  our  non-dimensionalization  the 
advection-diffusion  equation  in  (1).  In  particular,  we  mapped  u  to  the  non-dimensional  fluid 
velocity  field  u  t— >  tt/||n||  and  £  to  the  non-dimensional  molecular  diffusivity  e  »  e/(^||it||), 
where  i  is  the  maximum  cell  period  and  ||tt||  =  (|,u|2)1//2  is  the  Hilbert  space  norm  of  u.  For 
the  fluid  velocity  fields  in  equation  (99),  V  =  [0,  27 r]2  and  l  =  2tt.  When  u  in  equation  (99) 
is  non-random  then  the  underlying  Hilbert  space  is  =  L2(V,m),  where  m  denotes  the 
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normalized  Lebesgue  measure  (uniform  distribution)  on  restricted  to  V,  and  ||  •  ||  denotes 
the  -norm,  with 


2  B2  +  C2 

\UBC  =  - X - 


\\uce\\2 


1  +  A2 
2 


(100) 


In  Section  III  B  we  discussed  the  effective  parameter  problem  for  the  setting  of  randomly 
perturbed,  periodic  flows.  For  the  sake  of  brevity,  we  consider  here  only  the  randomly 
perturbed  cat’s  eye  flow,  with  the  parameter  A  uniformly  distributed  on  the  interval  [0, p] , 
having  second  moment  p2  / 3.  Numerically,  it  is  natural  to  set  =  L2(m  x  P)  [26],  where 
P  is  the  probability  measure  associated  with  the  random  variable  A.  In  this  case,  by 
Fubini’s  theorem  [18],  (•)  can  be  considered  to  denote  spatial  averaging  followed  by  statistical 
averaging  and  the  formula  for  ||wce||2  in  (100)  is  replaced  by 


\UCE\ 


3  +  p2 

6 


(101) 


We  now  discuss  in  more  detail  our  discrete,  matrix  formulation  of  the  effective  parameter 
problem.  To  illustrate  how  to  generalize  these  ideas  to  dimension  d  larger  than  d  =  2,  we 
will  maintain  aspects  of  our  general  notation.  In  this  discrete  setting,  the  spatial  region 

V  =  [0,  27r]d,  for  example,  is  replaced  by  a  square  lattice  Vd  of  size  L  containing  Ld  equally 
spaced  points  in  V.  As  discussed  in  the  beginning  of  Section  IV  A,  the  differential  operators 

V  and  V-  are  replaced  by  finite  difference,  matrix  operators  V  and  —  Vr,  respectively,  with 
suitable  boundary  conditions.  Periodic  boundary  conditions  will  be  assumed  throughout 
this  section.  Since  these  matrices  operate  on  vectors,  the  d-dimensional  lattice  Vd  must  be 
bijectively  mapped  to  a  one  dimensional  lattice  Vn  of  size  N  =  Ldd.  The  specific  structure 
of  Vjv  and  V  depend  on  the  bijective  mapping  0  chosen.  In  our  computations  discussed  in 
this  section,  we  used  the  mapping  0  described  in  [37]. 

Under  the  mapping  0,  a  spatially  dependent  d-dimensional  vector  held  v(x),  say,  is 
bijectively  mapped  to  a  discretized  constant  vector  of  length  TV  that  contains  all  of  the 
spatial  information  of  v(x)  for  x  E  Vd.  For  example  the  spatially  dependent  vector  v(x)  = 
(vi(x) ,  V2(x))  is  mapped  to  (ui,u2),  where  the  vectors  v^,  i  =  1,2,  are  constant  and  of 
length  Ld.  Similarly,  the  d-dimensional  standard  basis  vector  e±  =  (1,0,...,  0)  is  mapped 
to  the  TV-dimensional  vector  (1,0, ...  ,0),  where  1  and  0  are  vectors  of  ones  and  zeros  of 
length  Ld,  respectively,  and  similarly  for  the  ej  for  j  —  2, . . . ,  d.  Therefore,  the  vectors  e3, 
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j  —  1, ...  ,d,  satisfying 

ej  =  0(ej)/Ld/2,  ej  ■  ek  =  Sjk,  (102) 

serve  as  “lattice  standard  basis  vectors.”  In  previous  sections,  we  deferred  the  description 
of  these  vectors  to  the  present  section  and,  for  simplicity,  used  the  notation  e3.  Now  that 
the  specific  nature  of  these  vectors  has  been  discussed,  we  will  henceforth  use  the  notation 
in  equation  (102).  With  this  convention,  the  division  by  Ld /2  in  (102)  takes  care  of  the 
L-scaling  in  discrete  representations  of  spatial  integrals,  where  the  normalized  Lebesgue 
measure  da?  becomes  the  discrete  differential  da?  =  (2n / L)d / {2n)d  when  V  =  [0,  27r]rf  and 
the  spatial  average  (£(a?)  £(a?))y  becomes  £>-C,/Ld.  As  another  example,  under  the  bijective 
mapping,  the  2x2  matrix  H  in  (97)  becomes  a  N  x  N  antisymmetric  banded  matrix,  where 
the  stream  function  'I' (a?)  is  represented  by  a  diagonal  Ld  x  Ld  matrix  and  the  zero  element 
0  is  represented  by  a  matrix  of  zeros.  In  higher  dimensions  d>  2  the  discrete  representation 
of  the  matrix  H  is  also  banded.  As  in  previous  sections,  for  notational  simplicity,  we  will  not 
make  a  notational  distinction  for  the  matrix  H  between  the  continuum  and  discrete  settings. 

In  Section  VIII  we  demonstrated  that,  for  the  discrete  setting,  the  spectral  measure  jijk 
underlying  the  Stieltjes  integral  representation  of  S*k  is  given  by 

d/b'fc(A)  =  ^  (mjk(n)5x i(dA)),  (103) 

n :  A*  <A 

where  A*,  n  —  1, . . . ,  Ad,  are  the  eigenvalues  of  the  matrix  —  ?U(  HUi  in  (94),  while  vari¬ 
ous  equivalent  representations  of  the  spectral  weights  rrijk(n),  j,k  =  1 , ,d,  are  given  in 
equation  (96).  In  our  computations  of  fijk,  we  used 

mjk(n)  =  Re  [  ([rjjtllf  H%)  (K]tUf  Hefc)  ],  n  =  1, . . . ,  Ku  (104) 

which  follows  from  equations  (38),  (95),  and  (96).  We  stress  that,  for  notational  simplicity, 
in  this  section  we  denote  p,jk  =  Re  Hjk,  j,  k  —  1,2,  as  indicated  in  equation  (104).  In  (104), 
r*,  n  =  1,  are  the  complex  eigenvectors  of  the  matrix  —  ?.UfHUi.  Consequently, 

pkk  is  a  positive  measure  and  pjk  is  a  signed  measure  for  j  ^  k.  The  size  of  the  matrix 
V  is  N  x  Ld,  where  N  =  Ldd.  For  d  =  2  the  nullity  of  V  is  1,  therefore,  the  size  of  Ui  is 
N  x  ( Ld  —  1),  so  that  the  Hermitian  matrix  — 1\)\  HUi  is  of  size  K\  =  Ld  —  1. 

We  denote  the  spectral  weights  rrijk(n)  associated  with  the  decomposition  fijk  =  n+k  —  pjk 
in  (30)  by  rn+k (n)  and  m~k(n),  where  m±k(n)  >  0.  Moreover,  we  also  denote  [18]  by  the 
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functions  [S^2]+  and  [S*2 


(Si2]  +  (£)  =  max{SJ2(e),0},  [S*12]  (e)  =  max{-S*2(e),  0},  (105) 

for  each  0  <  £  <  oo,  so  that  S^2  =  [S*2]  +  —  [S*2]“,  [S^2]±(e)  =  S *2(e;  /if2),  and  [S^2]  =*=  >  0. 

In  the  case  of  a  non-random  fluid  velocity  held  u,  we  used  L  =  200  so  that  K\  =  39,  999. 
The  eigenvalues  A4  and  eigenvectors  r\  of  the  non-random  Hermitian  matrix  —  ?UfHUi 
were  computed  using  the  Matlab  function  eig().  In  this  case,  the  averaging  (•)  in  (103) 
is  interpreted  as  spatial  averaging  over  the  period  cell  V.  In  the  random  setting,  we  used 
L  =  100  so  that  K\  =  9,  999.  In  this  case,  the  averaging  (•)  in  (103)  is  interpreted  as  spatial 
averaging  followed  by  ensemble  averaging  over  ~  103  statistical  trials. 

The  numerical  accuracy  of  the  eigenvalue  problem  is  determined  by  the  eigenvalue  con¬ 
dition  numbers  /C(A4 ),  n  —  1, . . . ,  Ad,  which  are  the  reciprocals  of  the  cosines  of  the  angles 
between  the  left  and  right  eigenvectors.  Large  eigenvalue  condition  numbers  of  a  Hermitian 
matrix  implies  that  it  is  near  a  matrix  with  multiple  eigenvalues,  while  eigenvalue  condition 
numbers  close  to  1  imply  that  the  eigenvalue  problem  is  well-conditioned.  The  eigenvalue 
problem  for  the  matrix  —  zllf  HUi  is  extremely  well  conditioned  with  max,,  1 1— /C(A4 )  |  ~  10“14 
typical,  which  were  computed  using  Matlab’s  function  condeig(). 

To  our  knowledge,  Matlab  does  not  provide  a  function  that  describes  the  accuracy  of 
the  computed  SVD  of  the  matrix  V  =  UEV7  .  In  order  to  better  understand  the  numerical 
accuracy  in  the  entries  of  the  matrix  U,  which  is  central  to  our  computational  method,  we 
performed  the  following  tests.  For  the  case  of  Dirichlet  boundary  conditions,  the  matrix 
V  is  of  full  rank,  hence  the  matrix  Laplacian  V7  V  is  invertible.  We  computed  the  matrix 
T  =  V(V7  V)-1VT  directly  using  Matlab’s  mldivide  function  A\B,  i.e.,  T  =  V((V7  V)\VT), 
and  also  using  the  SVD  of  the  matrix  V,  with  T  =  UU7  .  We  then  computed  the  component¬ 
wise  maximum  difference  rnax/m  [V((V7  V)\VT)  —  UUT]/mj.  When  L  =  100  and  L  =  200 
this  difference  is  ~  10~15,  which  gives  an  idea  of  the  accuracy  of  the  SVD  of  V  for  the 
rank  deficient,  periodic  setting.  In  all  of  our  computations,  Matlab’s  sparse  architecture  was 
employed  wherever  possible  to  reduce  roundoff  error. 
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B.  Numerical  Results 

We  now  discuss  our  numerical  results.  It  was  shown  [3]  in  the  continuum  setting  that  for 
shear  flow  in  the  x-direction,  the  measure  gn  is  given  by  a  5-measure  concentrated  at  the 
spectral  origin,  while  fi 22  =  0,  and  similarly  for  shear  flow  in  the  //-direction.  As  a  baseline 
result,  we  computed  the  spectral  measures  and  effective  diffusivities  for  BC-shear-flow  in 
both  the  x  and  //-directions,  which  are  obtained  for  parameter  values  ( B ,  C )  =  (0, 1)  and 
(. B,C )  =  (1,0),  respectively.  Our  computations  of  the  components  jijk,  j,k  =  1,2,  of  the 
spectral  measure  for  BC-shear-flow  displayed  in  Fig.  2  are  in  agreement  with  the  theoretical 
prediction  in  [3]. 

Fig.  2  displays  the  streamlines  for  BC-shear-flow  in  (a)  the  x-direction  and  (b)  the  y- 
direction.  In  Fig.  2c  the  components  S*fc,  j,k  =  1,2,  of  the  effective  diffusivity  tensor  are 
displayed  for  BC-shear-flow  in  the  x-direction.  The  analogous  result  for  BC-shear-flow  in 
the  //-direction  is  visually  identical  to  Fig.  2c  under  the  mapping  S*x  S22,  he.,  under 

the  exchange  of  the  colors  black  blue.  The  components  y,jk,  j,k  =  1,2,  of  the  spectral 
measure  are  displayed  for  BC-shear-flow  in  (d)  the  x-direction  and  (e)  the  //-direction. 

We  focus  our  discussion  on  the  results  for  BC-shear-flow  in  the  x-direction,  as  the  dis¬ 
cussion  regarding  BC-shear-flow  in  the  //-direction  is  analogous.  For  all  n  —  1, ...  ,K±,  the 
spectral  weights  m 22(71)  in  Fig.  2d  associated  with  the  //-direction  satisfy  77122(71)  <  1CT29, 
while  771  ^"2  (n)  <  lO-28  in  the  bulk  of  the  spectrum  with  a  peak  near  the  spectral  origin  with 
spectral  weights  satisfying  rn^2  (n)  <  1CT16.  With  the  effects  of  finite  resolution  L  <  00  as 
well  as  numerical  errors  in  the  computed  components  of  the  matrix  — illf  HUi  and  its  eigen¬ 
value  decomposition,  associated  with  roundoff  error  due  to  a  machine  epsilon  of  ~  10“16, 
these  spectral  weights  can  be  considered  “numerically  zero.”  The  spectral  weights  for  the 
x-direction  satisfy  77111(77)  <  10-28  in  the  bulk  of  the  spectrum,  while  the  weights  near  the 
spectral  origin  satisfy  1CT9  <  77in  <  1CT1,  resembling  a  5-measure  with  virtually  all  of  its 
mass  concentrated  near  the  origin.  This  is  consistent  with  theoretical  predictions  [3].  Due 
to  the  antisymmetry  of  the  real-valued  matrix  Iff  HUi,  its  complex  eigenvectors  and  purely 
imaginary  eigenvalues  come  in  complex  conjugate  pairs  [22],  Consequently,  the  eigenvalues 
of  the  Hermitian  matrix  — zllf  HUi  come  in  positive-negative  pairs  with  identical  spectral 
weights,  resulting  in  the  symmetry  about  the  y- axis  displayed  by  the  spectral  measures 
in  Fig.  2. 
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FIG.  2.  Shear  flow  baseline  result.  The  streamlines  of  BC-shear-flow  in  (a)  the  x-direction  and  (b) 
the  y-direction.  (c)  The  e  behavior  of  our  computations  of  the  components  S*k,  j,k  =  1,  2,  of  the 
effective  diffusivity  corresponding  to  shear  flow  in  the  x-direction.  The  spectral  weights  rrijk  of  the 
spectral  measure  Re  fijk  for  shear  flow  in  (d)  the  x-direction  and  (e)  the  y-direction.  Consistent 
with  theoretical  predictions,  the  measure  associated  with  the  direction  of  the  flow  resembles  a  delta 
measure  centered  at  the  origin,  while  the  other  two  components  have  spectral  weights  rrijk  with 
very  small  magnitudes. 


Due  to  the  high  concentration  of  measure  mass  in  fin  very  near  the  spectral  origin,  our 
computation  of  S*x  displayed  in  Fig.  2c  behaves  like  it’s  being  governed  by  a  delta  function 
concentrated  at  the  origin.  In  particular,  Fig.  2c  shows  that  the  computed  £  behavior  of 
displayed  in  black  color  with  solid  line-style  lays  right  on  top  of  the  graph  of  its  upper  bound 
e  [1  +  /U°fc/e2]  given  in  (29),  with  ~  2.5204  x  10~2,  displayed  in  black  color  and  dash-dot 
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line-style.  (We  had  to  increase  the  line-width  of  the  upper  bound  to  be  able  to  distinguish 
between  the  two  black  lines.)  Due  to  the  extremely  small  magnitudes  of  the  spectral  weights 
m2 2  and  mf2%  with  measure  masses  /i22  ~  2.7050  x  10-30,  [/x?2]+  ~  5.2119  x  10~18,  and 
[yi2]~  ~  1-6910  x  10-16,  the  upper  and  lower  bounds  for  S22  and  S*2  hi  equations  (29) 
and  (31)  are  very  close  to  £  and  0,  respectively;  The  graph  of  S22  is  virtually  right  on  top  of 
the  lower  bound  £  in  cyan  color  and  solid  line-style,  and  the  magnitudes  of  [S*2]  +  and  [S^2] — 
are  so  small  they  don’t  even  appear  on  the  graph.  Since  the  support  of  the  spectral  measure 
is  contained  in  the  interval  [—1, 1],  the  components  of  the  effective  diffusivity  approach  their 
bare  molecular  diffusivity  value  eSjk  for  large  e. 

In  [38]  we  developed  Fourier  methods  for  the  rigorous  computation  of  the  spectral  measure 
fijk  for  BC- cell  flow,  with  B  —  C  —  1.  In  particular,  the  eigenvalue  problem  Aipn  =  A nipn  as¬ 
sociated  with  the  operator  A  =  A_1[w*V]  was  transformed  into  an  infinite  algebraic  system 
of  equations,  defining  a  discrete,  generalized  eigenvalue  problem.  The  Fourier  coefficients  of 
the  eigenfunctions  (pn,  n  —  1,  2,  3, . . .,  for  the  continuum  setting  comprised  the  components 
of  the  generalized  eigenvectors  in  the  discrete  setting.  Moreover,  motivated  by  the  theoret¬ 
ical  findings  in  the  current  work,  we  provided  a  rigorous  extension  of  the  results  given  here 
to  the  setting  of  a  time-dependent  fluid  velocity  field,  where  A  =  A-1[cA  +  tt«V]  and  dt 
denotes  partial  differentiation  in  time.  Furthermore,  we  used  abstract  methods  of  functional 
analysis  to  generalize  Lemma  1  to  the  continuum,  steady  and  dynamic  settings.  The  Fourier 
methods  in  [38]  were  also  generalized  to  the  setting  of  a  time-dependent  fluid  velocity  field. 
Since  we  already  treated  BC- cell  flow  in  [38],  and  for  the  sake  of  brevity,  we  now  turn  our 
attention  to  a  discussion  regarding  our  results  for  “cat’s  eye  flow.” 

Since  the  streamlines  for  cat’s  eye  flow  in  Fig.  1  are  symmetric  about  the  line  y  =  x,  as 
discussed  above,  we  anticipate  that  yu  =  /x22.  Our  computations  of  the  components  /i^, 
j,  k  =  1,2,  of  the  spectral  measure  displayed  in  Figures  3  and  4,  for  non-random  A,  indicate 
that  this  is  indeed  the  case.  A  closer  look  at  these  figures  reveals  a  deeper  symmetry,  namely 
that  /in  =  M22  =  |/ii2 1 ,  where  |/xi2 1  =  yf2  +  A12  is  the  total  variation  of  the  signed  measure 
yi2  introduced  in  equation  (30). 

Since  the  operator  A  =  A_1[V‘tt]  is  compact  [8],  its  spectra  is  discrete  with  a  limit  point 
at  the  spectral  origin  A  =  0  [49].  This  limit  point  behavior  of  the  measures  fijk ,  j,  k  =  1,2, 
can  be  seen  in  all  of  the  panels  of  Fig.  3.  When  the  parameter  A  =  0,  the  streamlines  of 
cat’s  eye  flow  are  closed  cell  structures,  as  shown  in  Fig.  1,  so  that  large  scale  transport 
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FIG.  3.  Transition  away  from  cat’s  eye  cell  flow.  The  spectral  weights  nijk  for  the  components 
Re  fijk,  j,k  =  1,  2,  of  the  spectral  measure  are  displayed  with  increasing  values  of  the  free  parameter 
A  from  left  to  right.  As  the  parameter  A  increases,  the  streamlines  of  the  flow  transition  away 
from  cell  structures  to  open  channels.  This  is  reflected  in  the  measure  by  a  dramatic  increase  in 
the  magnitude  of  the  spectral  weights  rrijk  associated  with  the  limit  point  of  the  measure  at  A  =  0, 
while  the  other  weights  change  only  slightly. 

occurs  only  when  e  >  0  [15].  In  this  case,  the  magnitude  of  the  spectral  weights  rrikk(n) 
and  m±k(n),  n  —  1, . . . ,  Ki,  associated  with  this  limit  point  at  A  =  0  are  <  10  28 ,  as  shown 
in  Fig.  3.  When  A  >  0,  open  channels  connect  neighboring  cells  and  large  scale  transport 
takes  place  both  in  thin  boundary  layers  and  within  the  channels  [15].  This  is  reflected  in 
the  spectral  measure  by  a  dramatic  increase  in  the  magnitude  of  the  spectral  weights  mfcfc(n) 
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FIG.  4.  Transition  toward  cat’s  eye  shear  flow.  The  spectral  weights  rrijk  for  the  components 
Re  fijk ,  i,  k  =  1,2,  of  the  spectral  measure  are  displayed  with  increasing  values  of  the  free  parameter 
A  from  left  to  right.  As  the  value  of  the  parameter  A  increases,  the  streamlines  become  more 
elongated  in  the  x-y  diagonal  direction,  becoming  shear  flow  when  A  =  1.  This  is  reflected  in 
the  spectral  measure  by  an  increase  in  the  breadth  of  the  spectral  region  with  significant  measure 
mass. 

and  m±k(n)  associated  with  the  limit  point  at  A  =  0,  by  more  than  15  orders  of  magnitude 
for  a  change  of  only  10-6  in  the  magnitude  of  A ,  while  the  spectral  weights  associated  with 
the  bulk  of  the  spectrum  change  only  subtlety,  as  shown  in  Fig.  3. 

As  the  value  of  A  increases  into  the  range  (10— 2 , 10°),  the  limit  point  near  A  =  0  persists. 
However,  the  increase  in  the  magnitudes  of  the  associated  spectral  weights  stops,  with  values 
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in  the  range  (lCrn,lCr5)  for  all  A  G  (1CT2,100).  In  this  regime,  the  limit  point  of  the 
spectral  measure  is  closely  surrounded  by  spectrum  in  the  bulk  having  spectral  weights  with 
comparable  magnitudes.  As  A  increases  in  the  range  (10_1,10°),  a  significant  transitional 
behavior  arises  in  the  spectral  measure  in  the  bulk  of  the  spectrum,  as  shown  in  Fig.  4.  In 
particular,  a  plateau  forms  in  the  spectral  measure  for  A  G  (—1,  —0.5)  U  (0.5, 1)  with  spectral 
weights  having  magnitudes  <  10-13.  Another  feature  also  arises  that  has  a  more  significant 
influence  on  the  behavior  of  the  effective  diffusivity  S*.  Namely,  the  appearance  of  spectra 
A  G  (—0.5, —0.2)  U  (0.2,  0.5)  with  measure  masses  in  the  range  (1CF10, 1CT2),  as  shown  in 
Fig.  4.  This  broadening  of  the  region  having  spectral  weights  with  magnitudes  as  large  as 
1CT2  from  A  G  (—0.2,  0,2),  present  for  all  A  G  [0,1],  to  A  G  (—0.5,  0,5)  has  a  significant 
influence  on  the  behavior  of  the  effective  diffusivity  S*,  as  shown  in  Fig.  5. 

Our  computations  of  the  components  S*fc,  j,  k  =  1,  2,  for  cat’s  eye  flow  are  displayed  in 
Fig.  5  along  with  their  upper  bounds  given  in  the  same  color  and  dash-dot  line-style.  Since 
the  support  of  fijk  is  contained  in  the  interval  [—1,1],  the  components  S*k  of  the  effective 
diffusivity  approach  their  bare  molecular  diffusivity  value  sSjk  for  large  e.  We  discussed 
above  that  our  computations  of  the  components  / ijk ,  j  =  1,  2,  of  the  spectral  measure  display 
the  symmetry  /in  =  /i22  =  | A*i2 1  -  Since  the  behavior  of  the  /ijk  govern  the  behavior  of  the 
corresponding  components  of  the  effective  diffusivity  S*fc,  the  symmetry  /in  =  /i22  =  I//12 1 
between  the  measures  gives  rise  to  the  symmetry  Sn(e)  =  Si^e)  =  £  +  X(e]  n\2)  +A"(e;  /x]~2) 
between  the  components  of  the  effective  diffusivity,  where  we  have  denoted  by  A"(e;z/)  = 
f  dz/(A)/(e2  +  A2),  e.g.,  S^(e)  =  e  +  X(e;  /in).  The  symmetry  =  S22  can  be  clearly  seen 
in  our  computations  of  S*fc,  j,  k  =  1,  2,  displayed  in  Fig.  5;  The  two  curves  lay  right  on  top 
of  one  another,  as  do  their  upper  bounds  as  /in  =  /i22.  The  lower  bounds  for  S£fc,  k  =  1,  2, 
have  been  omitted  in  the  figure  panels,  as  they  are  virtually  right  on  top  of  the  looser  lower 
bound  £  (displayed  in  cyan  color)  due  to  the  small  measure  masses  <  10-3.  We  have  also 
numerically  explored  the  approximate  relationship  Sn  ~  £  +  [S^]"1"  +  [Si2]~  by  plotting  S*n 
S^,  and  £  +  [S*12]+  +  [S*2]_  on  one  graph.  For  most  values  of  A  and  £  considered,  the  three 
curves  lay  virtually  on  top  of  each  other  (not  shown),  and  when  there  is  a  deviation  of 
£  +  [Si2]+  +  [S;j-  from  Sn,  it  is  slight. 

Recall,  we  demonstrated  in  Fig.  3  that  for  A  G  (0, 1CT2),  the  limit  point  A  =  0  of  the 
spectral  measure  Hjk  has  weights  rrijk  with  magnitudes  that  increase  dramatically  as  A 
increases  from  zero,  with  magnitudes  in  the  interval  (1CT11, 10”5)  when  A  ~  1CT2.  However, 
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FIG.  5.  Transitional  behavior  of  the  effective  diffusivity  from  cat’s  eye  cell  flow  to  shear  flow.  The 
behavior  of  the  components  S*k,  j,k  =  1,  2,  of  the  effective  diffusivity  as  a  function  of  the  molecular 
diffusivity  e  and  increasing  values  of  the  parameter  A  from  left  to  right  and  top  to  bottom.  The 
upper  bounds  corresponding  to  S*k  are  in  dash-dot  line-style  and  are  the  same  color  for  Skk  and 
red  for  S*2,  while  the  lower  bound  £  for  Skk  is  in  cyan  color  and  solid  line-style.  For  small  values  of 
the  parameter  A,  the  enhancement  in  the  effective  diffusivity  is  more  pronounced  for  small  values 
of  £.  However,  as  the  value  of  the  parameter  A  increases  and  the  flow  transitions  from  cell  to  shear 
structure,  there  is  a  substantial  enhancement  in  the  effective  diffusivity  for  larger  values  of  e.  This 
is  due  to  the  behavior  of  the  spectral  measure  discussed  in  Figures  3  and  4. 

in  the  bulk  of  the  spectrum,  the  magnitudes  of  the  spectral  weights  change  only  subtlety. 
Consequently,  for  A  G  (0, 1CT2)  this  transitional  behavior  of  the  spectral  measure  fikk  governs 
primarily  the  small  £  behavior  of  S*kk,  as  shown  in  the  panels  of  Fig.  5  corresponding  to  A  =  0 
and  A  =  10~2.  The  transitional  behavior  of  S^2  is  more  pronounced  due  to  the  lack  of  the 
eSjk  term  for  j  ^  k. 

Recall,  we  demonstrated  in  Fig.  4  that  when  the  parameter  A  increases  in  the  interval 
(lO-1, 10°),  spectra  A  G  (—0.5, —0.2)  U  (0.2,  0.5)  appear  with  measure  masses  in  the  range 
(10~10, 10'2).  This  broadens  the  influence  of  the  spectral  measure  fijk  over  the  effective 
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FIG.  6.  Migration  of  positive  measure  mass  from  the  limit  point  at  the  spectral  origin.  As  the 
free  parameter  A  of  cat’s  eye  flow  increases  from  zero,  the  magnitude  of  the  spectral  weights 
increase  dramatically.  Moreover,  the  spectra  associated  with  the  positive  weights  mf2  migrate 
away  from  the  spectral  origin  until  the  limit  point  is  comprised  only  of  negative  valued  mass. 

diffusivity  S*k  to  larger  values  of  e,  greatly  enhancing  it  above  its  bare  molecular  diffusivity 
value  e5jk,  as  shown  in  the  panels  of  Fig.  5  for  A  >  0.1.  Note  that,  since  gn  =  p 22  =  | M12 1 , 
we  have  the  inequalities  S*x  >  £  +  [SJ2]+  and  S^,  >  e  +  [S];2]_,  with  S^,  =  S22,  as  shown  in 
Fig.  5. 

When  the  figure  panels  associated  with  // 1 2  in  Fig.  3  are  plotted  in  log-log  scale  as  shown 
in  Fig.  6,  the  following  is  revealed.  As  A  increases  from  zero,  the  spectra  of  the  limit  point 
with  positive  measure  mass  migrates  away  from  A  =  0,  so  that  the  limit  point  eventually 
consists  of  spectra  with  only  negative  measure  mass.  Consequently,  as  e  decreases  below 
10~3,  this  influence  of  p\2  on  S*2  becomes  more  dominant  and  S*2  changes  sign,  becoming 
negative,  as  shown  in  Fig.  5.  As  the  value  of  £  approaches  the  location  of  this  limit  point, 
the  numerical  approximation  breaks  down  due  to  effects  of  finite  resolution  L. 

We  now  discuss  our  computations  of  the  components  p,jk  and  S*k,  j,k  =  1,2,  of  the 
spectral  measure  and  effective  diffusivity,  respectively,  for  cat’s  eye  flow  with  random  pa¬ 
rameter  A  uniformly  distributed  on  the  interval  [0,  p] .  For  each  statistical  trial  of  a  sample 
space  of  ~  103  statistical  trials  and  a  system  resolution  L  =  100,  we  computed  ev¬ 
ery  eigenvalue  A3  and  eigenvector  r3,  n  =  1, . . . ,  Ad,  of  the  matrix  — 1\}\  HUi  to  form  the 
spectral  measure  pjk  in  equation  (103).  In  order  to  visually  determine  the  behavior  of  the 
function  Pjk(^)  =  (Q(A)ej,  e*,)  underlying  the  spectral  measure  pjk,  we  plot  a  histogram 
representation  of  p,jk{ A)  called  the  spectral  function,  which  we  will  also  denote  by  pjk( A). 
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FIG.  7.  Spectral  functions  and  effective  diffusivities  for  randomly  perturbed  cat’s  eye  flow.  The 
random  parameter  A  is  uniformly  distributed  on  the  interval  [0. p] .  The  spectral  functions  Pjk{ A) 
are  displayed  with  corresponding  effective  diffusivities  S*fc  directly  below  for  various  values  of  p. 
increasing  from  left  to  right.  As  p  increases  and  the  streamlines  of  the  flow  become  more  elongated 
in  the  x-y  direction,  on  average,  the  region  about  the  spectral  origin  A  =  0  with  substantial 
measure  mass  increases  in  breadth  and  magnitude.  This  gives  rise  to  a  substantial  enhancement  in 
the  components  S*k  of  the  effective  diffusivity  for  larger  values  of  the  molecular  diffusivity  e.  The 
color  scheme  of  the  panels  for  S*fc  is  the  same  as  that  in  Fig.  5. 


We  now  describe  how  we  computed  this  graphical  representation  of  the  measure  Hjk ■  First, 
the  spectral  interval  /  A  E  was  divided  into  V  sub-intervals  Iv,  v  —  1, . . . ,  V,  of  equal  length 
1/V.  Second,  for  fixed  v,  we  identified  all  of  the  eigenvalues  that  satisfy  A^(u;)  G  Iv,  for 
n  —  1, . . . ,  K\  and  uj  e  Oq-  The  assigned  value  of  Hjk( A)  at  the  midpoint  A  of  the  interval 
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Iv,  is  the  sum  of  the  spectral  weights  rrijk(ou)  associated  with  all  such  A h(u>)  G  Iv.  In  our 
computations  of  the  spectral  functions,  we  typically  used  V  ~  102. 

Consistent  with  the  symmetries  of  the  random  flow,  our  computations  of  the  spectral 
function  satisfies  pn(A)  =  /422(A),  hence  the  ensemble  averaged  components  S*k  of  the 
effective  diffusivity  also  satisfy  S*x  =  S22,  as  shown  in  Fig.  7.  Similar  to  onr  computations 
for  non-random  A,  when  p  =  0.1  the  measure  mass  of  j,  k  —  1,2,  near  the  spectral  origin 
is  quite  small  and,  on  average,  the  region  with  significant  magnitude  increases  in  breadth  as 
p  increases.  This  average  increase  in  the  breadth  of  the  region  with  significant  mass  gives 
rise  to  a  substantial  enhancement  of  the  components  S*fc  of  the  effective  diffusivity  above 
the  bare  molecular  diffusivity  values  sSjk-  Our  results  for  cat’s  eye  flow  with  random  A  also 
demonstrate  the  influence  of  resolution  L  in  the  numerical  computations.  In  particular,  with 
a  decrease  in  resolution  L  from  L  =  200  in  Figures  3  and  4  to  L  —  100  in  Fig.  7,  we  see  that 
the  accuracy  of  the  numerical  computations  break  down  for  £  10  3  with  L  =  100  instead 

of  £  ~  10“4  with  L  =  200,  indicated  by  a  1/e  divergence.  For  the  continuum  setting,  the 
limit  point  of  the  spectrum  at  A  =  0  can  be  discrete  with  finite  or  infinite  multiplicity,  and 
can  even  be  continuous  [51]. 


X.  CONCLUSIONS 

In  Section  III  we  adapted  and  extended  a  method  [2,  3]  that  provides  the  rigorous  Stielt- 
jes  integral  representations  for  the  symmetric  S*  and  antisymmetric  A*  parts  of  effective 
diffusivity  tensor  D*  shown  in  equation  (28).  These  integral  representations  involve  the 
molecular  diffusivity  £  and  a  spectral  measure  [ijk  of  a  self-adjoint  operator  that  acts  on  the 
Hilbert  space  of  curl-free  vector  holds.  In  Section  IV  we  considered  a  discrete  formulation 
of  Section  III  for  the  case  that  the  matrix  Laplacian  is  of  full  rank,  hence  invertible.  A 
matrix  analysis  showed  that  the  spectral  measure  is  given  in  terms  of  the  eigenvalues  and 
eigenvectors  of  a  Hermitian  matrix.  In  Section  IV  B  we  provided  a  projection  method  which 
revealed  that  many  of  the  spectral  weights  of  the  discrete  spectral  measure  are  identically 
zero,  while  the  others  are  determined  by  a  much  smaller  Hermitian  matrix.  This  method 
stabilizes  and  increases  the  efficiency  of  numerical  computations  of  pjk,  which  enables  more 
accurate  computations  of  the  effective  diffusivity  tensor  D*. 

In  Section  V  we  returned  to  the  continuum  setting,  adapting  and  extending  a  different 


53 


method  [8,  41]  that  provides  the  rigorous  Stieltjes  integral  representations  in  equation  (28), 
involving  the  molecular  diffusivity  £  and  a  (possibly  different)  spectral  measure  /.q^  of  a  self- 
adjoint  operator  that  acts  on  a  Sobolev  space  of  scalar  fields.  In  Section  VI  we  considered 
a  discrete  formulation  of  Section  V  for  the  case  that  the  matrix  Laplacian  is  of  full  rank. 
A  matrix  analysis  showed  that  the  spectral  measure  is  given  in  terms  of  the  generalized 
eigenvalues  and  eigenvectors  associated  with  a  pair  of  Hcrmitian  matrices  of  the  same  size 
as  that  arising  in  the  above  mentioned  projection  method. 

In  Section  VII  we  used  properties  of  the  singular  value  decomposition  of  the  matrix  gra¬ 
dient  V  to  reveal  symmetries  between  the  two  discrete  approaches  formulated  in  Sections  IV 
and  VI,  establishing  in  Lemma  1  that  the  two  approaches  yield  equivalent  spectral  repre¬ 
sentations  of  the  effective  diffusivity  tensor  D*  when  the  matrix  Laplacian  is  of  full  rank. 
In  particular,  we  established  in  the  proof  of  Lemma  1  that  the  eigenvalues  and  generalized 
eigenvalues  underlying  the  spectral  measures  for  each  method  are  in  fact  eigenvalues  of  a 
Hermitian  matrix  arising  in  both  methods.  Moreover,  the  eigenvectors  wn  and  generalized 
eigenvectors  zn  of  the  two  methods  are  related  by  wn  =  Vzn,  which  leads  to  the  equivalence 
of  the  discrete  spectral  measures  of  the  two  approaches. 

In  Section  VIII  we  generalized  Lemma  1  to  the  case  that  matrix  Laplacian  is  rank  defi¬ 
cient,  hence  non-invert ible.  This  extends  the  numerical  methods  developed  in  Sections  IV 
and  VI  to  the  setting  of  periodic  boundary  conditions,  for  which  the  matrix  Laplacian  is 
singular.  In  over  25  years  since  the  formulation  [2,  3]  of  Stieltjes  integral  representations 
for  the  effective  diffusivity  tensor  D*,  analytical  calculations  of  D*  have  been  obtained  for 
only  a  few  simple  flows,  such  as  shear  flow.  Our  results  of  Section  VIII  help  overcome  this 
limitation  by  providing  a  mathematical  foundation  for  rigorous  computation  of  D*. 

In  Section  IX  we  employed  the  numerical  methods  formulated  in  Section  VIII  to  com¬ 
pute  the  components  Sjk,  j,  k  —  1, . . . ,  d,  of  S*  for  some  model  2D  periodic  flows,  by  directly 
computing  the  associated  spectral  measure  Re  /.q^.  As  a  baseline  result,  we  computed  Sjk 
and  Re  /q*,  for  B C-shear-flow,  for  which  the  spectral  measure  is  known  [3].  Our  numerical 
results  are  in  good  agreement  with  the  theoretical  result.  We  also  computed  the  transitional 
behavior  of  Sjk  and  Re  /q*,  for  “cat’s  eye”  flow  as  a  function  of  a  free  parameter,  both  for 
the  random  and  non-random  settings.  Consistent  with  the  symmetries  of  the  flow,  our  com¬ 
putations  indicate  that  Re  qn  =  Re  /i22-  Our  computations  of  Re  /q^,  j,k  =  1,2,  for  cat’s 
eye  flow  also  reveal  a  deeper  symmetry,  namely  Re  fi\\  =  Re  q22  =  |Re  /q2|,  where  |Re  /q2| 
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is  the  total  variation  of  the  measure  Re  fi  12.  Our  computations  of  S jj.  are  consistent  with 
these  symmetries,  as  well  as  rigorous  bounds  that  we  derived  from  the  analytic  properties 
of  the  Stieltjes  integral  representation  for  S jk  and  the  associated  measures. 

Motivated  by  the  theoretical  findings  in  the  current  work,  in  [38]  we  provided  a  rigorous 
extension  of  the  results  given  here  to  the  setting  of  a  time-dependent  fluid  velocity  field 
u  =  u(t,x).  Furthermore,  we  used  abstract  methods  of  functional  analysis  to  generalize 
Lemma  1  to  the  continuum,  steady  and  dynamic  settings.  We  are  currently  exploring  the 
extension  of  these  methods  to  the  time-stochastic  setting,  which  is  relevant  to  atmospheric 
and  oceanic  flows. 
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Appendix  A:  Properties  of  the  linear  operator  A 

In  this  section  we  derive  various  properties  of  the  linear  operator  A  =  A_1[tt-V]  defined 
in  equation  (48).  In  particular,  we  demonstrate  that  A  is  antisymmetric  on  the  Hilbert 
space  7 i1-2  defined  in  (46).  Moreover,  we  show  that  A  is  bounded  on  H1:2  and  we  provide 
an  upper  bound  for  || ^4 1| i;2  when  u  is  uniformly  bounded  on  the  period  cell  V. 

We  first  show  that  the  incompressibility  condition  V-w  =  0  implies  that  the  operator 
A  is  antisymmetric  on  H1'2  [8],  (Af,h) i)2  =  —(f,Ah)i:2-  On  the  Hilbert  space  H  defined 
in  (45)  the  linear  operator  A-1  satisfies  (AA-1/,  h)  =  (/,  h)  in  a  distributional  sense,  for  all 
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f,h(zH  [17,  33].  Consequently,  for  all  /,  h  e  7i1,2  we  have 

W,/»)1>2  =  ([V(A-1)(u-V)/]-V/r)  (Al) 

=  -<[(«•  V)/],*) 

=  -([V.(u/)],/i) 

=  (/,[(«.  V^]) 

=  (/,[A(A-1)(u-V)/i]) 

=  -(V/.[V(A-1)(u-V)/r]) 

=  -(f,Ah)  il2. 


Now,  we  derive  the  bound  on  ||A||  given  in  equation  (52).  From  the  Cauchy-Schwartz 
inequality  |(/,  h)\  <  ||/||  ||/i||  we  have  that 

WAfWh  =  |(V[A-1(u.V/)]-V[A-1(u.V/)])|  (A2) 

=  \-([A-\u.Vf)},(u.Vf))\ 

<  ||A_1('u-V/)||  ||t*.V/|| 

<  llA_1ll  ll«-v/||2. 


We  now  provide  an  upper  bound  for  ||if  V/||  when  the  components  Uk,  k  =  1, . . . ,  d,  of  the 
fluid  velocity  field  u  are  uniformly  bounded  on  the  period  cell  V.  By  the  Cauchy-Schwartz 
inequality,  |£»£|  <  |£|  |£|,  we  have 


\\u.Vf\\2  =  (\u-Vf\2)  <  <|u|2  |  V/|2)  <  sup|w(a;)|2  ||/||22. 

x£V 


The  result  in  equation  (52)  is  now  clear. 


(A3) 


Appendix  B:  Derivation  of  equation  (87) 

In  this  section  we  provide  a  derivation  of  equation  (87).  Equation  (85)  allows  equa¬ 
tion  (80)  to  be  written  as  [Vi£iRi](di  +  zAi^ViEiRi]^-  =  Uj.  This  and  equation  (83) 
imply  that 


(Bl) 


VfXj  =  Z^Riieh  +  zAO^RlE^Vf^. 
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This  formula  and  V  =  IRZRVf  imply  that 

VXj  =  UiRi(eli  +  (B2) 

Therefore,  since  U^IR  =  li  and  R{Ri  =  R,  we  clearly  have  the  first  formula  in  (87)  for  S*k 
with  Z\  =  RjEj^Vf.  From  equation  (85)  we  have  that 

(VTHVX/Xt)  =  <»[ViS1R1]A1[RiE1Vf]5c3..Xfc)  =  (tp.RARlEJVf^-Vfx,,}.  (B3) 

This  formula,  R|Ri  =  R,  Sf  =  Si,  and  equation  (Bl)  clearly  imply  the  second  formula  for 
A*fc  in  equation  (87),  with  Z\  =  R}SR1Vf. 
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